我用的是 gdal.Translate
通过掩码提取光栅。两种光栅数据共享相同的CRS和相同的像素大小。下面是一段代码。
def by_mask(in_dir, in_fname, mask_dir, mask_fname, out_dir, out_fname):
mask_ds = gdal.Open(os.path.join(mask_dir, mask_fname), GA_ReadOnly)
mask_proj = mask_ds.GetProjection()
mask_geotrans = mask_ds.GetGeoTransform()
ulx = mask_geotrans[0]
uly = mask_geotrans[3]
lrx = ulx + mask_geotrans[1] * mask_ds.RasterXSize
lry = uly + mask_geotrans[5] * mask_ds.RasterYSize
ds = gdal.Open(os.path.join(in_dir, in_fname), GA_ReadOnly)
gdal.Translate(destName=os.path.join(out_dir, out_fname),
srcDS=ds,
options=gdal.TranslateOptions(
format='GTiff',
projWin=[ulx, uly, lrx, lry],
options=['-co', 'COMPRESS=DEFLATE'],
xRes=res_info(in_dir=mask_dir, in_fname=mask_fname)[0],
yRes=res_info(in_dir=mask_dir, in_fname=mask_fname)[1]
))
我得到了错误信息 ERROR 1: Error: Computed -srcwin 11577 6424 1236 -873 has negative width and/or height.
有什么建议吗?
这是我的解决方案。欢迎提出任何建议或问题。
from osgeo import gdal
from gdalconst import GA_ReadOnly
import os
import osr
def by_mask(in_dir, in_fname, mask_dir, mask_fname, out_dir, out_fname,
mask_rev_lst):
"""
extract raster data with another mask raster data, working well with
inconsistent extents
:param in_dir: directory of target raster data
:type in_dir: str
:param in_fname: filename of target raster data
:type in_fname: str
:param mask_dir: directory of mask raster data
:type mask_dir: str
:param mask_fname: filename of mask raster data
:type mask_fname: str
:param out_dir: directory of raster data to be saved
:type out_dir: str
:param out_fname: filename of raster data to be saved
:type out_fname: str
:param mask_rev_lst: values to masked
:type mask_rev_lst: list
:return:
:rtype:
"""
src_f = os.path.join(in_dir, in_fname)
mask_f = os.path.join(mask_dir, mask_fname)
dest_f = os.path.join(out_dir, out_fname)
src_ds = gdal.Open(src_f)
src_nodata = src_ds.GetRasterBand(1).GetNoDataValue()
exists_nodata = False if src_nodata is None else True
mask_ds = gdal.Open(mask_f)
mask_arr = mask_ds.ReadAsArray()
mask_geotrans = mask_ds.GetGeoTransform()
mask_proj = mask_ds.GetProjection()
mem_driver = gdal.GetDriverByName('Mem')
mem_ds = mem_driver.Create('', mask_ds.RasterXSize, mask_ds.RasterYSize,
1, gdal.GDT_Float32)
mem_ds.SetProjection(mask_proj)
mem_ds.SetGeoTransform(mask_geotrans)
mem_srs = osr.SpatialReference()
mem_srs.ImportFromWkt(mask_proj)
if src_nodata is not None:
mem_ds.GetRasterBand(1).SetNoDataValue(src_nodata)
gdal.Warp(destNameOrDestDS=mem_ds,
srcDSOrSrcDSTab=src_ds,
options=gdal.WarpOptions(
resampleAlg=gdal.gdalconst.GRA_NearestNeighbour,
srcNodata=src_nodata,
dstNodata=src_nodata
))
mem_arr = mem_ds.ReadAsArray()
del mem_ds
mem_arr_lst = []
for el in mask_rev_lst:
temp_arr = mem_arr.copy()
temp_arr[~(mask_arr == el)] = 1e29
if exists_nodata:
mem_arr_lst.append(temp_arr)
else:
mem_arr_lst.append(temp_arr[~(temp_arr == 1e29)])
dest_driver = gdal.GetDriverByName('GTiff')
dest_ds = dest_driver.Create(dest_f,
mask_ds.RasterXSize, mask_ds.RasterYSize,
1, gdal.GDT_Float32,
options=['COMPRESS=DEFLATE'])
dest_ds.SetGeoTransform(mask_geotrans)
dest_ds.SetProjection(mask_proj)
dest_ds.GetRasterBand(1).WriteArray(*mem_arr_lst)
if src_nodata is not None:
dest_ds.GetRasterBand(1).SetNoDataValue(1e29)
dest_ds.FlushCache()