1.重合两个不同年份的tif的分辨率import rasterio import numpy as np from rasterio.warp import reproject, Resampling from rasterio.windows import from_bounds base_file rG:\研究区ldmod2020rectangleNew.tif hx_file rG:\ldmod2020hx.tif output_file rG:\ldmod2020rectangleNEWhx.tif # # 1. 读取底图 # with rasterio.open(base_file) as base: base_data base.read(1) profile base.profile.copy() base_crs base.crs base_transform base.transform base_width base.width base_height base.height base_bounds base.bounds print(底图:) print( size:, base_width, base_height) print( CRS:, base_crs) print( resolution:, base.res) print( bounds:, base_bounds) # # 2. 读取河西走廊土地利用 # with rasterio.open(hx_file) as hx: hx_data hx.read(1) hx_crs hx.crs hx_transform hx.transform hx_bounds hx.bounds print(\n河西走廊:) print( size:, hx.width, hx.height) print( CRS:, hx_crs) print( resolution:, hx.res) print( bounds:, hx_bounds) # # 3. 检查范围 # if ( hx_bounds.left base_bounds.left or hx_bounds.right base_bounds.right or hx_bounds.bottom base_bounds.bottom or hx_bounds.top base_bounds.top ): raise ValueError( ldmod2020hx.tif 并没有完全包含在 ldmod2020rectangle.tif 范围内。 ) # # 4. 计算 HX 在底图中的位置 # window from_bounds( hx_bounds.left, hx_bounds.bottom, hx_bounds.right, hx_bounds.top, transformbase_transform ) window window.round_offsets().round_lengths() row_start int(window.row_off) row_end row_start int(window.height) col_start int(window.col_off) col_end col_start int(window.width) print(\n替换范围:) print( rows:, row_start, row_end) print( cols:, col_start, col_end) # # 5. 创建与底图窗口完全一致的 HX 数据 # target_height row_end - row_start target_width col_end - col_start target_transform rasterio.windows.transform( window, base_transform ) hx_aligned np.full( (target_height, target_width), 255, dtypenp.uint8 ) # # 6. 重投影 最近邻重采样 # reproject( sourcehx_data, destinationhx_aligned, src_transformhx_transform, src_crshx_crs, dst_transformtarget_transform, dst_crsbase_crs, resamplingResampling.nearest, src_nodata255, dst_nodata255 ) # # 7. 替换 # base_data[ row_start:row_end, col_start:col_end ] hx_aligned # # 8. 输出 # profile.update( driverGTiff, heightbase_height, widthbase_width, transformbase_transform, crsbase_crs, dtypeuint8, count1, nodata255, compresslzw ) with rasterio.open(output_file, w, **profile) as dst: dst.write(base_data, 1) print(\n) print(融合完成) print() print(输出:) print(output_file) print(\n输出数据:) print(size:, base_width, base_height) print(CRS:, base_crs) print(resolution:, profile[transform].a, abs(profile[transform].e)) print(unique:, np.unique(base_data))2.融合两个tif为新的tif数据import rasterio import numpy as np from rasterio.warp import reproject, Resampling from rasterio.windows import from_bounds # -------------------------------------------------- # 输入 # -------------------------------------------------- base_file rG:\ldmod2020rectangle_Resample1.tif hx_file rG:\ldmod2020hx.tif # 输出 out_file rG:\ldmod2020rectangleNew.tif # # 2. 打开底图和 HX 数据 # with rasterio.open(base_file) as base, \ rasterio.open(hx_file) as hx: print( BASE ) print(CRS:, base.crs) print(Size:, base.width, base.height) print(Resolution:, base.res) print(Bounds:, base.bounds) print(NoData:, base.nodata) print(\n HX ) print(CRS:, hx.crs) print(Size:, hx.width, hx.height) print(Resolution:, hx.res) print(Bounds:, hx.bounds) print(NoData:, hx.nodata) # # 3. 检查坐标系 # if base.crs ! hx.crs: raise ValueError( f两个 TIFF 坐标系不同\n fBASE: {base.crs}\n fHX: {hx.crs} ) # # 4. 检查 HX 是否位于底图范围内 # if hx.bounds.left base.bounds.left: raise ValueError(HX 左边界超出底图范围) if hx.bounds.right base.bounds.right: raise ValueError(HX 右边界超出底图范围) if hx.bounds.bottom base.bounds.bottom: raise ValueError(HX 下边界超出底图范围) if hx.bounds.top base.bounds.top: raise ValueError(HX 上边界超出底图范围) # # 5. 计算 HX 在 BASE 中对应的窗口 # window from_bounds( hx.bounds.left, hx.bounds.bottom, hx.bounds.right, hx.bounds.top, transformbase.transform ) # 四舍五入到完整像元 window window.round_offsets().round_lengths() print(\n 替换窗口 ) print(window) # # 6. 读取 BASE # base_data base.read(1) # # 7. 读取 HX # hx_data hx.read(1) # # 8. 创建与 BASE 对齐的 HX 数组 # target_height int(window.height) target_width int(window.width) hx_aligned np.full( (target_height, target_width), 255, dtypenp.uint8 ) # HX 在 BASE 中对应区域的 transform target_transform rasterio.windows.transform( window, base.transform ) # # 9. 最近邻重采样 # reproject( sourcehx_data, destinationhx_aligned, src_transformhx.transform, src_crshx.crs, dst_transformtarget_transform, dst_crsbase.crs, src_nodata255, dst_nodata255, resamplingResampling.nearest ) # # 10. 只替换 HX 有效值 # row_start int(window.row_off) row_end row_start target_height col_start int(window.col_off) col_end col_start target_width base_region base_data[ row_start:row_end, col_start:col_end ] # HX 有效像元 valid hx_aligned ! 255 print(\nHX有效像元数量:, np.sum(valid)) print(HX总像元数量:, hx_aligned.size) # 只用 HX 有效值替换 base_region[valid] hx_aligned[valid] # 放回底图 base_data[ row_start:row_end, col_start:col_end ] base_region # # 11. 保存结果 # profile base.profile.copy() profile.update( dtyperasterio.uint8, count1, nodata255, compresslzw ) with rasterio.open(out_file, w, **profile) as dst: dst.write(base_data, 1) print(\n) print(处理完成) print(输出文件) print(out_file) print()3.将新的tif转化为能被WPS识别的二进制数据。#!/usr/bin/env python3 import rasterio import numpy as np from pathlib import Path # # 1. 输入 TIFF # tif_file ldmod2020hx.tif # # 2. 输出 WPS GEOG 数据目录 # out_dir Path(landuse2020) out_dir.mkdir(parentsTrue, exist_okTrue) # # 3. 读取 GeoTIFF # with rasterio.open(tif_file) as src: data src.read(1) width src.width height src.height transform src.transform crs src.crs nodata src.nodata print() print(Input TIFF information) print() print(width :, width) print(height :, height) print(dtype :, data.dtype) print(CRS :, crs) print(NoData :, nodata) print(transform :, transform) print(unique values:) print(np.unique(data)) # -------------------------------------------------------- # 检查数据 # -------------------------------------------------------- if width ! 1085 or height ! 576: print(WARNING: raster size is different from expected.) if crs.to_epsg() ! 4326: raise ValueError(ERROR: TIFF is not EPSG:4326) # # 检查土地利用类别 # valid data[data ! nodata] print(\nValid value range:) print(min , valid.min()) print(max , valid.max()) print(\nValid classes:) print(np.unique(valid)) # -------------------------------------------------------- # 土地利用数据使用 uint8 # -------------------------------------------------------- data data.astype(np.uint8) # # WPS geogrid要求 # 从南向北写入 # # GeoTIFF: # 第0行 北边 # # WPS: # 第1行 南边 # # 所以需要上下翻转 # data_wps np.flipud(data) # # 4. 写 WPS binary # binary_file out_dir / 00001-01085.00001-00576 data_wps.tofile(binary_file) print(\nBinary file created:) print(binary_file) # # 5. 计算 WPS 网格参数 # dx transform.a dy abs(transform.e) # GeoTIFF左上角是外边界 west transform.c north transform.f # 第一个像元中心 lon_first west dx / 2 # 最南边第一个像元中心 lat_first north - (height - 0.5) * dy print(\n) print(WPS parameters) print() print(dx , dx) print(dy , dy) print(known_lon , lon_first) print(known_lat , lat_first) # # 6. 创建 index 文件 # index_file out_dir / index index_content ftypecategorical category_min1 category_max21 projectionregular_ll dx{dx:.15f} dy{dy:.15f} known_x1 known_y1 known_lat{lat_first:.15f} known_lon{lon_first:.15f} wordsize1 tile_x{width} tile_y{height} tile_z1 missing_value255 unitscategory descriptionMODIS IGBP 21-category land use 2020 Hexi Corridor mminluMODIFIED_IGBP_MODIS_NOAH iswater17 islake21 isice15 isurban13 with open(index_file, w) as f: f.write(index_content) print(\nIndex file created:) print(index_file) print(\n) print(DONE) print()