gdal reprojectimage 截断图像

Question

我正在使用 gdal python 重新投影 grib 文件。我这样做 python 而不是 gdalwarp 是因为我不想承担编写新文件并在 python 脚本中读回以进行进一步处理的打击。运行 gdalwarp 非常慢并且会创建非常大的文件。我目前正在使用 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 )

问题是它似乎切断了我的数据。我不确定如何更改重新投影图像的范围以包含所有数据。这是它的样子： result

它看起来像什么 expected

Answer 1

我必须找到 max/min X 和 max/min Y。用它代替 lr、ul 点。

编辑：基本上你不能在计算新的地理变换时使用右下角、左上角的点，因为当 reprojected/transformed 在你的角边界（lr，ul）之外时，你的网格中间可能有一个点。因此，您需要重新投影沿边界的所有点并获得最小值和最大值 (x,y)。然后使用这些最小值和最大值来计算新的地理变换。

gdal reprojectimage 截断图像

gdal reprojectimage cuts off image

python

gdal