2

I have developed a semantic segmentation method to map certain objects in aerial imagery. Throughout this project I download pictures (4800x4800 pxs) where I know the min_x, min_y, max_x and max_y in world coordinates, e.g. (6.212454957892032, 51.58908266914109, 6.219403242107967, 51.59339941284951). Resolution is 10 cm per pixel.

To then segment the pictures I crop the large picture (4800x4800) into multiple pictures with a resolution of 320x320. To identify each of those cropped pictures I calculated its center coordinate as follows (from upper left to right then to south):

dlat = (side * 360) / (2 * np.pi * r) #dlat for 32m in degrees
side = 32
r = 6371000 #avg earth radius 

minx = float(minx)
miny = float(miny)
maxx = float(maxx)
maxy = float(maxy)

identifier = (minx, miny, maxx, maxy)

# Takes a 4800x4800 image tile and returns a list of 320x320 pixel images

tile = np.array(tile)

images = []
coords = []

N = 0
S = 4800
W = 0
E = 4800

# y coordinate is dlat/2 degrees, i.e. 16 meters, south of the maximum y coordinate.
y_coord = maxy - dlat/2

while N < S:

    W = 0

    x_coord = minx + (((side * 360) / (2 * np.pi * r * np.cos(np.deg2rad(y_coord))))/2) #16 m to the middle

    while W < E:
        # The first image is taken from the upper left corner, we then slide from left
        # to right and from top to bottom

        images.append(tile[N:N + 320, W:W + 320])
        coords.append((x_coord, y_coord))

        x_coord += (((side * 360) / (2 * np.pi * r * np.cos(np.deg2rad(y_coord)))))

        W = W + 320

    N = N + 320

    y_coord = y_coord - dlat

As I also use this approach to then convert the identified pixel polygons to real world coordinates I run into minor differences between the underlying map and the identified objects (10-20m difference) sometimes minor. So I ask myself what is wrong with my approach and where could I make it more efficient by using proper packages.

I actually did some research but it is rather confusing to me because I do not use GeoTIFF files or similar. I get a numpy array as output which I then rasterize with rasterio to retrieve the polygons in pixel coordinates. I convert these using the following function to real world coordinates (using the middle point of each picture as an identifier and fix point to calculate the other coordinates):

def centroid_coord(center_coord, distance_vector, size = 320):
'''
returns lat, lon of array centroid
:param center_coord: tuple of center coords, identifier of picture tile
:param distance_vector: tuple of distance vector, measured from center point in px
:param size: scaling factor for distance from center point in pixels (half the size of the segmented output (320)
:return: lat,lon array centroid
'''
dist_px_x, dist_px_y = distance_vector
x_center, y_center = center_coord
x_min = x_center - (((32 * 360) / (2 * np.pi * r * np.cos(np.deg2rad(y_center))))/2)
y_new = y_center + (side/size) * dist_px_y * (dlat/side)
x_new = x_min + (side/size) * ((size/2) + dist_px_x) * 360 * (1/(2 * np.pi * r * np.cos(np.deg2rad(y_new))))
return (x_new, y_new)

I think that this approach is not optimal and that there is a much more correct and efficient way using proper packages.

How can I do this?

Array: Binary array with a size (320,320) indicating where a certain object is or not. Using rasterio finding polygons of the objects in pixel coordinates. Output: Convert polygons in pixel coordinates to world coordinates using the identifier of the input image (center point) My problems:

  • Is the center point calculation (world-coordinatewise) correct like this?
  • Should I take the upper left corner point of the image to ease the process?
  • How can I convert the polygons with pixel coordinates correctly to world coordinate polygons?
2

Although the length of the arc of the upper parallel and the length of the arc of the lower parallel of the limits of the original tile, measured on a sphere of 6371 km radius, measure 480 m with a variation less than 10 centimeters that a pixel measures, In general we should consider it a mistake to establish a spatial resolution in linear dimensions and think of everything else as angular.

Let's say better that the spatial resolution is 0.10m x 0.10m the pixel, and that the coordinate system is flat. In 480 m around, I am a terraplanista. I measure with precision instruments to 480 m around and draw a plane of the terrain, not only without considering the terrestrial curvature, but also without considering the reduction to the ellipsoid (not to mention the geoid undulation or the deviation of the vertical).

Let's say the system is flat and tangent to the sphere in the center of the image. Coordinates from the center can be considered as eastings and northings in meters of a topocentric CRS. Or let's say the system is a cylinder, which axis rests in the plane of the equator and is tangent to the spheric surface at the center of the image. Both options are valid in 480 m around, calculable with pyproj, and as easy to use with the sphere as with the WGS84 ellipsoid. How to do it depends of your pyproj version, which can differ if your are using a conda environment. In any case, how to do it would be another answer.

About your code, and trying to leaving aside the wrong assumption that each pixel can have the same spatial resolution if defined by geographic coordinates:

Here: y_new = y_center + (side/size) * dist_px_y * (dlat/side), when you have an even pixel matrix, the center is in the corner of four of them. Then, the distance to the first pixel to any side is the size of that pixel / 2, the distance to other pixels can be generalized as : (coordinate of the pixel - 0.5) * size of pixels.

And here: x_new = x_min + (side/size) * ((size/2) + dist_px_x) * 360 * (1/(2 * np.pi * r * np.cos(np.deg2rad(y_new)))), I don't know why you go to x origin with the length of the center parallel arc but then to the right with the length of the current parallel arc. Each pixel up will have a difference of longitude.

I think that nothing of that may produce a 16 m error, wich is like the distance from the center to the edge of a submatrix. And the problem may be in the distance_vector calculation.

About all other formulas, I have not checked them but understand your intention and agree with the underlying math used to accomplish it.

4
  • Thanks for the explanation! I will look into pyproj. I already checked it yesterday but I figured this cannot produce this kind of error. Regarding your points: No thats fine in python. It definitely is a 320x320. What do you mean by size of pixels? The size of a pixel in meters? That is what I tried to do with side/size. How could I correct the calculation for x_new, e.g. with pyproj? distance_vector = [x - side/2, y - side/2] is how I calculate the distance vector. There might be a mistake in my thinking: should it be [side/2 - x, side/2 - y] that might explain the error. – inneb Mar 22 '20 at 9:36
  • You are welcome. It's ok, I hadn't use np.tile before. Yes, you have the size of pixel in meters, and then you multiply for the size of one meter in degrees, it's ok too, but there is a half of pixel that you must subtract to transform pixel coordinates to meters or degrees. With pyproj you define the CRS and just transform. You can do it from topocentric to geocentric and then to geographic. Or from a custom Transverse Mercator to geographic. Are x and y the pixel coordinate from the lower left corner in the submatrix? Consider: [x - side/2, -y + side/2] if upper left corner. – Gabriel De Luca Mar 22 '20 at 12:55
  • I am seeing it is not np.tile, it's a slicing and it's ok. I removed that part of my answer. Can you get pyproj 2.6.0 or 1.9.6? The idea is to calculate all submatrix centers and other relative coordinates in meters, with a 0.1m pixel size. In the last step return to the original source CRS, in meters centered in the center of the source tile, and transform from there to geographic coordinates. – Gabriel De Luca Mar 22 '20 at 13:28
  • x and y are the pixel coordinates of a point of the polygon. side/2 is my center point of the slice as I also use the center coordinate as an identificator of the slice. Such that I know how far the polygon point is away from the center of the slice.I then use the distance vector with the x_new / y_new transformation to obtain the coordinate of the polygon point in a world coordinates. Does it make any sense to you now? I get your thoughts on the transformation. I'll try that out. :-) Many thanks! – inneb Mar 22 '20 at 14:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.