参考https://www.jianshu.com/p/821c741ff169
拼接
def rasterwarp(base_file,inpath,outpath):
if not os.path.exists(outpath):
os.mkdir(outpath)
out_file = os.path.join(outpath,'mcd12q1.tif')
inputrasfile1 = gdal.Open(base_file, gdal.GA_ReadOnly) # 第一幅影像
inputProj1 = inputrasfile1.GetProjection()
options = gdal.WarpOptions(srcSRS=inputProj1, dstSRS=inputProj1, format='GTiff', resampleAlg=gdalconst.GRA_Max)
rasterlist = []
filelist = pathlib.Path(inpath).glob('*.tif')
for rasterfile in filelist:
rasterfile = str(rasterfile)
inputrasfile2 = gdal.Open(rasterfile, gdal.GA_ReadOnly) # 第二幅影像
rasterlist.append(inputrasfile2)
gdal.Warp(out_file, rasterlist, options=options)
裁剪
gdal.Warp('D:\MCD12Q1\globe\\arcpy\con_result\\final2_mcd12q1.tif', 'D:\MCD12Q1\globe\\arcpy\con_result\\re_mcd12q1.tif',outputBounds=[-180,-60,-160,60],dstSRS=Projection,resampleAlg=gdalconst.GRA_NearestNeighbour)
重采样
def warp_image(img_file: str, base_file: str, out_file: str) -> int:
src = gdal.Open(img_file, gdalconst.GA_ReadOnly)
src_proj = src.GetProjection()
# src_geotrans = src.GetGeoTransform()
nband = src.RasterCount
dataType = src.GetRasterBand(1).DataType
match_ds = gdal.Open(base_file, gdalconst.GA_ReadOnly)
match_proj = match_ds.GetProjection()
match_geotrans = match_ds.GetGeoTransform()
wide = match_ds.RasterXSize
high = match_ds.RasterYSize
dst = gdal.GetDriverByName('GTiff').Create(out_file, wide, high, nband, dataType)
dst.SetGeoTransform(match_geotrans)
dst.SetProjection(match_proj)
gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_NearestNeighbour)
del dst
return 0