gdal重投影切断图像

Question

我正在使用gdal python重新设计一个grib文件。我之所以做这个python而不是gdalwarp，是因为我不想受到编写新文件和阅读python脚本来进行进一步处理的影响。运行gdalwarp非常慢，并且会创建非常大的文件。我目前正在使用来自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 )

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

​

它看起来像什么？

​

Answer 1

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

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