outfile = r"D:\微信公众号\掩膜\匹配后数据.tif"
ref_ds = gdal.Open(ref_file)
base_xsize = ref_ds.RasterXSize
base_ysize = ref_ds.RasterYSize
base_proj = ref_ds.GetProjection()
base_gt = ref_ds.GetGeoTransform()
xmin = base_gt[0]
ymax = base_gt[3]
xmax = xmin + base_xsize * base_gt[1]
ymin = ymax + base_ysize * base_gt[5]
#待处理数据
pend_ds = gdal.Open(pend_file)
#数据匹配
dataset = gdal.Warp(outfile,
pend_ds,
width=base_xsize,
height=base_ysize,
srcNodata=-999 ,
dstNodata=-999,
dstSRS=str(base_proj),
outputBounds=(xmin, ymin, xmax, ymax),
format="GTiff",
resampleAlg=gdal.GRA_Bilinear)
import gdal
#匹配后文件
infile = r"D:\微信公众号\掩膜\匹配后数据.tif"
#掩膜文件
maskfile = r"D:\微信公众号\掩膜\基准数据.tif"
#结果文件
outfile = r"D:\微信公众号\掩膜\掩膜后数据.tif"
#读取匹配后文件
inds = gdal.Open(infile)
indata = inds.ReadAsArray()
cols = inds.RasterXSize
rows = inds.RasterYSize
geo = inds.GetGeoTransform()
proj = inds.GetProjection()
band = inds.GetRasterBand(1)
#读取掩膜文件
maskds = gdal.Open(maskfile)
maskdata = maskds.ReadAsArray()
band_mask = maskds.GetRasterBand(1)
mask_nodata = band_mask.GetNoDataValue()
indata[maskdata == mask_nodata] = -999
#输出结果，无效值设置为-999
driver = gdal.GetDriverByName("Gtiff")
out_ds = driver.Create(outfile, cols, rows, 1, band.DataType)
out_ds.SetGeoTransform(geo)
out_ds.SetProjection(proj)
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(indata)
out_band.SetNoDataValue(-999)
欢迎关注我的个人公众号：小Rserimport gdal#基准数据ref_file = r"D:\微信公众号\掩膜\基准数据.tif"#待匹配数据pend_file = r"D:\微信公众号\掩膜\待匹配数据.tif"# 结果文件outfile = r"D:\微信公众号\掩膜\匹配后数据.tif"ref_ds = gdal.Open(ref_file)base_xsize = ref_ds.RasterXSizebase_ysize = ref_ds.RasterYS.
如果是像元大小不统一造成的，我们只需要对这些栅格进行重采样即可，重采样使这些栅格像元大小
一致
后再进行
掩
膜
（裁剪）。（较容易）
如果是投影坐标系不统一造成的，就需要在
掩
膜
（裁剪）之前对其重新投影栅格使投影坐标系统一。
着重讲下第二个解决办法：投影坐标系不统一
使用如下矢量
数据
作为
掩
膜
（裁剪）的范围
文章目录背景解决方案
不同
数据
源做重采样或者投影变换后，用同一
掩
膜
进行提取，往往会出现范围不
一致
的情况。这里说的范围不
一致
是指图像的上下左右平面坐标不完全
一致
，会出现细微的差别，导致栅格像素无法完全重叠。本文教大家如何解决该问题。
我调研过rasterio，发现并没有解决方案。calculate_default_transform函数无法固定输出图像范围，只能固定输入图像范围，而输入图像因为不同源，往往是没法固定的。
我的解决方案是引入arcpy包。
arcpy是arcmap中的toolbo
现有栅格
数据
以及矢量
数据
（红色框部分）
本来想用矢量
数据
对栅格
数据
进行
掩
膜
提取，得到研究区栅格
数据
，过程如下，在ArcGIS中AecToolbox>Spatical Analyst Tools>Extraction>Extract by Mask,
过程很简单，运行也正常，但结果却出现错误，
数据
...
