Pyrhon ＋ GDAL

Raster を複数使用した安定計算。

解像度の違いでエラーを吐きます。
Arc で ポリゴンからラスターへ変換していたのですが、解像度が完全には一致しないようでした。仕方がないので、リサンプリングすることに。
以前、GDAL を使って Raster を読み込んでいたので、そこに追記。

# Open Raster
def open_src(file,band):
    src = gdal.Open(file, gdal.GA_ReadOnly)
    if src is None:
        print('file not found')
    else:
        print(src.GetDescription() )
        xy=src.GetGeoTransform()#
        print(xy) 
        print(src.GetProjection())
        print(src.GetMetadata() )
        arr_b1=src.GetRasterBand(band).ReadAsArray(0)
        print(arr_b1.shape)
        nodata=src.GetRasterBand(band).GetNoDataValue())
        print(nodata)
        return src,xy,arr_b1,nodata


#図化
def plot_raster(arr_b1, nodata):
    arr_b1 = np.where(arr_b1==nodata, np.nan, arr_b1)
    print('max:',np.nanmax(arr_b1))
    print('min:',np.nanmin(arr_b1))
    print(arr_b1.shape)    
    plt.imshow(arr_b1)
    plt.colorbar()
    return arr_b1

     
#測地系変換
#EPSG：4612 から 4326に変換して保存
ds = gdal.Warp(file_out, file_in, srcSRS="EPSG:4612", dstSRS="EPSG:4326")
ds=None


#2つのラスターの領域、解像度を（ほぼ）揃える
#切り出し領域を決める
minx=136
miny=35
maxx=137
maxy=36
#領域 をforium で確認
fmap = folium.Map(location=[(miny+maxy)/2,(minx+maxx)/2], zoom_start=8)
line = [[miny,minx],[miny,maxx],[maxy,maxx],[maxy,minx],[miny,minx]]
folium.PolyLine(locations=line).add_to(fmap)
fmap
#ピクセル数を決める
width=18000
height=11100
#resampling
cut_file0="./output/cut0.tif" 
cut_file1="./output/cut1.tif" 
ds=gdal.Translate(cut_file0, file_in0, projWin=[minx,maxy,maxx,miny], width=width, height=height, noData=nodata)
ds=gdal.Translate(cut_file1, file_in1, projWin=[minx,maxy,maxx,miny], width=width, height=height, noData=nodata)
ds=None


#保存
def save_dst(rasterfile, arr_new, nodata, src):
    driver = gdal.GetDriverByName("GTiff")  
    dst = driver.Create(rasterfile, 
                       xsize=arr_new.shape[1],
                       ysize=arr_new.shape[0],
                       bands=1,
                       eType = gdal.GDT_Float32)
    dst.SetGeoTransform(src.GetGeoTransform())
    dst.SetProjection(src.GetProjection())
    band1=dst_ds.GetRasterBand(1)
    band1.WriteArray(arr_new)
    band1.SetNoDataValue(nodata)
    dst.FlushCache() 
 

Arc でラスターを2枚作成・書き出していましたので、マルチバンドにしてから1枚出力すれば解像度が完全に一致したのですが。
https://pro.arcgis.com/ja/pro-app/2.8/tool-reference/data-management/composite-bands.htm

ま、今後使えるかもしれませんので残しておきましょう。