Search code examples
pythongridgisgdal

Why does gdal_grid turn image upside-down?

I'm trying to use gdal_grid to make an elevation grid from a surface in a geojson. I use this command:

gdal_grid -a linear:radius=0 inputSurface.geojson outputFile.tif

It seems to give the correct pixel values, but if I open the result in Global Mapper or QGIS, the image is flipped/mirrored in a horizontal axis, such that the tif is directly below the surface and upside-down. What is the reason for this and how do I fix it??

Update

I already tried changing the geotransform, but it hasn't totally fixed my problem.

I looked at the resulting image in gdalinfo and found out that the upper left corner is actually the lower left corner, so I set it using the SetGeoTransform. This moved it to the correct location, but it is still upside-down. (This may by dependent on the projection, which might cause problems later)

I also tried looking at the pixel width in the geotransform as mentioned below:

Xgeo = GT[0] + Xpixel*GT[1] + Yline*GT[2]
Ygeo = GT[3] + Xpixel*GT[4] + Yline*GT[5]

The image returned by gdal_grid has a positive GT[5], but unfortunately changing it to -GT[5] doesn't change anything.

The code I used to change the geotransform:

transform = list(ds.GetGeoTransform())
transform = [upperLeftX, transform[1], 0, upperLeftY, 0, -transform[5]]
ds.SetGeoTransform(transform)
Solution

  • I ended up having more problems with gdal_grid, where it just crashes at seemingly random places, so I'm using the scipy.interpolate-function called griddata in stead. This uses a meshgrid to get the coordinates in the grid, and I had to tile it up because of memory restrictions of meshgrid.

    import scipy.interpolate as il #for griddata
import numpy as np
# meshgrid of coords in this tile
gridX, gridY = np.meshgrid(xi[c*tcols:(c+1)*tcols], yi[r*trows:(r+1)*trows][::-1])

## Creating the DEM in this tile
zi = il.griddata((coordsT[0], coordsT[1]), coordsT[2], (gridX, gridY),method='linear',fill_value = nodata) # fill_value to prevent NaN at polygon outline
raster.GetRasterBand(1).WriteArray(zi,c*tcols,nrows-r*trows-rtrows)

    The linear interpolation seems to do the same as gdal_grid is supposed to. This was actually effected by making the 5'th element in the geotransform negative as described in the question update.

    See description at scipy.interpolate.griddata.

    A few things to note:

    1. The point used in the geotransform should be upper-left
    2. The resolution in y-direction should be negative
    3. In the projection (at least the ones I use) positive y-direction is up
    4. In numpy arrays positive y-direction is down
    5. When using gdal's WriteArray it uses the upper left corner

    Hope this helps other people's confusion.