Pixel coordinate to world coordinates using Python

Question

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?

Answer 1

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.