Reproject coordinates with rasterio creates a irregular grid from an regular grid

Question

I have geotiff files load into xarray with a crs = EPSG:31467. I want to transform/reproject (don't know if there is a difference) these files into EPSG:4326. To do that, I use rasterio.warp.transform function which needs 1D arrays for x,y. To generate these i use numpy.meshgrid and flatten functions. Here is a small example with my data:

import numpy

#Longitude and Latitude in EPSG:31467
lon = [3280914, 3281914, 3282914]
lat = [6103001, 6102001, 6101001]

#create 2d meshgrid
xv, yv = np.meshgrid(lon, lat)


xv, yv 

    (array([[3280914, 3281914, 3282914],
        [3280914, 3281914, 3282914],
        [3280914, 3281914, 3282914]]),
 array([[6103001, 6103001, 6103001],
        [6102001, 6102001, 6102001],

        [6101001, 6101001, 6101001]]))

Now I have a sequence of different longitude [3280914, 3281914, 3282914] for the same latitude [6103001, 6103001, 6103001] When i now use rasterio.transform(src_crs, dst_crs, x, y) these sequences disappear and i dont unterstand why?!

from rasterio.warp import transform

# Compute the lon/lat coordinates with rasterio.warp.transform
lon, lat = transform('EPSG:31467','EPSG:4326',
                     xv.flatten(), yv.flatten())



np.asarray(lon).reshape(3,3), np.asarray(lat).reshape(3,3)

> (array([[5.57397386, 5.58957607, 5.6051787 ],
>         [5.57473921, 5.59033795, 5.60593711],
>         [5.57550412, 5.5910994 , 5.60669509]]),  array([[55.00756605, 55.00800488, 55.00844171],
>         [54.9985994 , 54.99903809, 54.99947477],
>         [54.98963274, 54.99007128, 54.99050782]]))

np.unique(xv).shape, np.unique(yv).shape

> ((3,), (3,))

np.unique(lon).shape, np.unique(lat).shape

> ((9,), (9,))

To change the reporjected coordinates back to xarray I have to get the same shape in sense of equality. Which process I don't understand, is it the function of transform or the concept of projections?

Answer 1

I can't understand what exactly you are trying to do after np.asarray(lon).reshape(3,3)

Which process I don't understand, is it the function of transform or the concept of projections?

It seems like you don't understand both.

EPSG:31467 and EPSG:4326 are fundamentally different types of data. EPSG:31467 is actually a planar rectangular coordinate system in zonal projection. EPSG:4326 is not a projection at all, it is a pure geodetic coordinates in WGS-84 terrestrial coordinate system with WGS-84 ellypsoid. What is exactly emportant here is that same coordinates in EPSG:31467 don't have to be same in EPSG:4326. Because in 4326 your coordinate is an angle and in 31467 your coordinate is a distance from equator or false meridien. Axes in these systems are not collinear and related with convergence of meridians parameter. So, if you change Norting or Easting in 31467, both latitude and logitute can change. Here you can notice an angle between blue lines (one cell is 31467 analogue) and black lines (whole grid is 4326 analogue) https://ru.wikipedia.org/wiki/%D0%A4%D0%B0%D0%B9%D0%BB:Soviet_topographic_map_kilometer_grid.svg

It's pretty easy to check, that transformation works correctly - just do it backwards.

lon, lat = transform('EPSG:31467','EPSG:4326',
                     xv.flatten(), yv.flatten())
x_check, y_check = transform('EPSG:4326', 'EPSG:31467', lon, lat)

#we'll have some troubles because of computational errors, so let's round
x_check = [int(round(i, 0)) for i in x_check]

print(lon)
print(x_check)
print(xv.flatten())

>[5.574033001416839, 5.5896346633743175, 5.605236748547687, 5.574797816145165, 5.5903960110246524, 5.605994628800234, 5.5755622060626155, 5.591156935778857, 5.6067520880717225]
>[3280914, 3281914, 3282914, 3280914, 3281914, 3282914, 3280914, 3281914, 3282914]
>[3280914 3281914 3282914 3280914 3281914 3282914 3280914 3281914 3282914]

Output examples that transform() returns you exactly what it expected to return.

Next code also works as it is expected (you can match output with above one):

print(np.asarray(lon).reshape(3,3))
print(xv)

>[[5.574033   5.58963466 5.60523675]
> [5.57479782 5.59039601 5.60599463]
> [5.57556221 5.59115694 5.60675209]]
>[[3280914 3281914 3282914]
> [3280914 3281914 3282914]
> [3280914 3281914 3282914]]

I have never worked with rasterio, so I can't provide you working solution. Some notes:

• I have no idea why do you need grid for raster transformation
• Rasterio docs are clear and have solution for you: https://rasterio.readthedocs.io/en/latest/topics/reproject.html#reprojecting-a-geotiff-dataset
• You can transform raster between crs directly. If not in rasterio, try osgeo.gdal (gdal.Warp(dst_file, src_file, srcSRS='EPSG:31467', dstSRS='EPSG:4326')
• Note the difference between reprojection and defining projection for raster. First changes image, second changes metadata. For correct work of direct transform, your GeoTIFF must have valid projection defenition in metadata (that matches actual projection of your raster)
• If you're not developing standalone app and just need to reproject 2-3 rasters, use QGIS and do it without coding. It's also helpfull to try understanding geodetic concepts on 2-3 examples in QGIS before coding. Just use it as a playground
• If you're not developing standalone app, you can solve your automatisation task in QGIS python API. You can test workflow with UI and then call some QGIS/GDAL tools from python script as batch. What is more - rasterio and all other packages will be avaluable for installation on QGIS' python. Of course, it's a bad idea for deployment unless you are creating a QGIS plugin
• In EPSG:31467 the coordinate value of 0.001 is 1 mm. So more precise is useless. In EPSG:4326 1 degree is 111.1 km approx (or 111.3*cos(lat)). So, you can calculate useful precise. Everything more than 4-5 digit after . may also be useless