读取由边界框定义的全局栅格的子集
打开覆盖地球的栅格并提取栅格的子集。
import gdal
# Path to a tiff file covering the globe
# http://visibleearth.nasa.gov/view.php?id=57752
tif_name = "/path_name/land_shallow_topo_21600.tif"
# Open raster in read only mode
ds = gdal.Open(tif_name, gdal.GA_ReadOnly)
# Get the first raster band
band = ds.GetRasterBand(1)
# Compute x/y resolution in degrees
resx = 360. / band.XSize
resy = 180. / band.YSize
# Define the geotransform used to convert x/y pixel to lon/lat degree
# [lon_topleft, lon_resolution, lat_skew, lat_topleft, lon_skew, lat_resolution]
geotransform = [-180, resx, 0.0, 90, 0.0, -1*resy]
# The inverse geotransform is used to convert lon/lat degrees to x/y pixel index
inv_geotransform = gdal.InvGeoTransform(geotransform)
# Define a longitude/latitude bounding box in degrees
# [lonmin, latmin, lonmax, latmax]
bbox = [-5, 40, 10, 55]
# Convert lon/lat degrees to x/y pixel for the dataset
_x0, _y0 = gdal.ApplyGeoTransform(inv_geotransform, bbox[0], bbox[1])
_x1, _y1 = gdal.ApplyGeoTransform(inv_geotransform, bbox[2], bbox[3])
x0, y0 = min(_x0, _x1), min(_y0, _y1)
x1, y1 = max(_x0, _x1), max(_y0, _y1)
# Get subset of the raster as a numpy array
data = band.ReadAsArray(int(x0), int(y0), int(x1-x0), int(y1-y0))
