Search code examples
pythongdal

gdal reprojectimage cuts off image


i'm using gdal python to reproject a grib file. i'm doing this python instead of gdalwarp is because i dont want to take the hit of writing a new file and reading back in the python script to do further processing. And running gdalwarp is very slow and creates very huge files. I'm currently using the reprojection code from http://jgomezdans.github.io/gdal_notes/reprojection.html .

osng = osr.SpatialReference ()
osng.ImportFromEPSG ( epsg_to )
wgs84 = osr.SpatialReference ()
wgs84.ImportFromEPSG ( epsg_from )
tx = osr.CoordinateTransformation ( wgs84, osng )
# Up to here, all  the projection have been defined, as well as a 
# transformation from the from to the  to :)
# We now open the dataset
g = gdal.Open ( dataset )
# Get the Geotransform vector
geo_t = g.GetGeoTransform ()
x_size = g.RasterXSize # Raster xsize
y_size = g.RasterYSize # Raster ysize
# Work out the boundaries of the new dataset in the target projection
(ulx, uly, ulz ) = tx.TransformPoint( geo_t[0], geo_t[3])
(lrx, lry, lrz ) = tx.TransformPoint( geo_t[0] + geo_t[1]*x_size, \
                                      geo_t[3] + geo_t[5]*y_size )
# See how using 27700 and WGS84 introduces a z-value!
# Now, we create an in-memory raster
mem_drv = gdal.GetDriverByName( 'MEM' )
# The size of the raster is given the new projection and pixel spacing
# Using the values we calculated above. Also, setting it to store one band
# and to use Float32 data type.
dest = mem_drv.Create('', int((lrx - ulx)/pixel_spacing), \
        int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)
# Calculate the new geotransform
new_geo = ( ulx, pixel_spacing, geo_t[2], \
            uly, geo_t[4], -pixel_spacing )
# Set the geotransform
dest.SetGeoTransform( new_geo )
dest.SetProjection ( osng.ExportToWkt() )
# Perform the projection/resampling 
res = gdal.ReprojectImage( g, dest, \
            wgs84.ExportToWkt(), osng.ExportToWkt(), \
            gdal.GRA_Bilinear )

The problem is that it seems to be cutting off my data. i'm not exactly sure how to change the extent of the reprojected image to contain the all of the data. here what it looks like: result

whats it suppose to look like expected


Solution

  • i had to find the max/min X and max/min Y. use that instead of the lr, ul points.

    Edit: Basically you cannot use the lower right, upper left points in calculating the the new geotransformation, since there could be a point in the middle of your grid that when reprojected/transformed is outside of your corner bounds (lr, ul). So you need to reproject all the points along your border and get the min and max (x,y). then use those min and max to calculate the new geotransform.