黄土高原有可能恢复成森林吗？沙漠会变成绿洲吗？

# 调用dde模型库中的Project Raster模型 import numpy as np from rasterio.warp import calculate_default_transform, reproject, Resampling from rasterio import crs import rasterio def pydde_rasterProject_run(src_file, EPSG_Code, tag_file): tag_crs = crs.CRS.from_epsg(EPSG_Code) with rasterio.open(src_file) as src_ds: profile = src_ds.profile tag_trans, tag_width, tag_height = calculate_default_transform( src_ds.crs, tag_crs, src_ds.width, src_ds.height, *src_ds.bounds) profile.update({ 'crs': tag_crs, 'transform': tag_trans, 'width': tag_width, 'height': tag_height, 'nodata': 0 }) with rasterio.open(tag_file, 'w', **profile) as tag_ds: for i in range(1, src_ds.count + 1): src_array = src_ds.read(i) tag_array = np.empty((tag_height, tag_width), dtype=profile['dtype']) reproject( source=src_array, src_crs=src_ds.crs, src_transform=src_ds.transform, destination=tag_array, dst_transform=tag_trans, dst_crs=tag_crs, resampling=Resampling.cubic, num_threads=2) tag_ds.write(tag_array, i) res = "success" return res input = input_dir + "N35E107_FABDEM_V1-0.tif" outDEM = output_dir + "DEM_P.tif" EPSG_Code = "4545" pydde_rasterProject_run(input, EPSG_Code, outDEM)

input = input_dir + "GBLU_v11_N35E107.shp" outGBLU = output_dir + "GBLU_P.shp" EPSG_Code = "EPSG:4545" # 调用gdal的ogr2ogr工具实现矢量数据重投影 command = 'ogr2ogr -f "ESRI Shapefile" -t_srs ' + EPSG_Code + ' ' + outGBLU + ' ' + input print(command) os.system(command)

outDEM = output_dir + "DEM_P.tif" outHS = temp_dir + "Hillshade.tif" wbt.hillshade( outDEM, outHS, azimuth=315.0, altitude=45.0, zfactor=None, callback=my_callback ) outDEM = output_dir + "DEM_P.tif" outHS = temp_dir + "Hillshade.tif" # 读取dem数据 dem = rxr.open_rasterio(outDEM) dem = dem[0] dem = dem.where(dem.values!=0) # Ignore nodata # 读取山地阴影图 hs = rxr.open_rasterio(outHS) hs = hs[0] hs = hs.where(hs.values!=0) # Ignore nodata alpha = 0.7 #将DEM叠加在山地阴影上，需要设置DEM图层透明度 fig, ax = plt.subplots(figsize=(12, 12)) hs.plot(cmap='gray', add_colorbar=False) # 先绘制灰色的山地阴影图作为底图 dem.plot(cmap='terrain', alpha=alpha, cbar_kwargs={'label': 'height'}) # 在山地阴影图上叠加半透明的dem # 隐藏xy坐标 ax.set_xticks([]) ax.set_yticks([]) # 添加标题 ax.set_title("DEM with hillshade") plt.show() # 清除变量，释放内存 del dem del hs plt.clf()

outDEM = output_dir + "DEM_P.tif" outFill = temp_dir + "Fill.tif" wbt.fill_depressions( outDEM, outFill, fix_flats=True, flat_increment=None, max_depth=None, callback=my_callback ) outFill = temp_dir + "Fill.tif" fill_image = rxr.open_rasterio(outFill) fill_image = fill_image[0] fill_image = fill_image.where(fill_image.values!=0) # Ignore nodata fig, ax = plt.subplots(figsize=(10, 10)) fill_image.plot(cmap='terrain', cbar_kwargs={'label': 'height'}) ax.set_xticks([]) ax.set_yticks([]) ax.set_title('Filled DEM data') plt.show() del fill_image plt.clf()

outFill = temp_dir + "Fill.tif" outFlowD = temp_dir + "FlowDir.tif" wbt.d8_pointer( outFill, outFlowD, esri_pntr=True, callback=my_callback ) outFlowD = temp_dir + "FlowDir.tif" fd_image = rxr.open_rasterio(outFlowD) fd_image = fd_image[0] fd_image = fd_image.where(fd_image.values!=-32768) # Ignore nodata fig, ax = plt.subplots(figsize=(10, 10)) # plt.imshow(fd_image,cmap='Spectral') # 对于取值范围为0-360°的流向数据，选用Spectral为色带 cmap = plt.cm.get_cmap('Spectral', 9) direct_cmap = ListedColormap([cmap(i) for i in range(9)]) boundaries = list([0, 1, 2, 4, 8, 16, 32, 64, 128]) norm = BoundaryNorm(boundaries, cmap.N, extend='max') fd_image.plot(cmap=direct_cmap, levels=boundaries, label='direction', cbar_kwargs={'label': 'direction'}) ax.set_xticks([]) ax.set_yticks([]) ax.set_title('Flow Direction') plt.show() del fd_image plt.clf()

outFlowD = temp_dir + "FlowDir.tif" outFlowAcc = temp_dir + "FlowAcc.tif" wbt.d8_flow_accumulation( outFlowD, outFlowAcc, out_type="cells", log=False, clip=False, pntr=True, esri_pntr=True, callback=my_callback ) outFlowAcc = temp_dir + "FlowAcc.tif" fa_image = rxr.open_rasterio(outFlowAcc) fa_image = fa_image[0] fa_image = fa_image.where(fa_image.values!=-32768) # Ignore nodata fig, ax = plt.subplots(figsize=(10, 10)) # 由于流量的绝对值相差过多，为了优化可视化效果，自定义对数坐标轴色带cbar norm = colors.LogNorm(vmin=fa_image.min(), vmax=fa_image.max()) im = ax.imshow(fa_image, norm=norm, cmap='Blues') cbar = plt.colorbar(im, ax=ax, format=ticker.ScalarFormatter(), ticks=[1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000], label='pixel number') cbar.ax.set_yticklabels(['1', '10', '100', '1000', '10000', '100000', '1000000', '10000000', '100000000']) # 设置刻度标签 ax.set_xticks([]) ax.set_yticks([]) ax.set_title('Flow Accumulation') plt.show() del fa_image plt.clf()

outFlowAcc = temp_dir + "FlowAcc.tif" threshold = '100000.0' outStream = temp_dir + "Stream.tif" wbt.extract_streams( outFlowAcc, outStream, threshold, zero_background=False, callback=my_callback ) outStream = temp_dir + "Stream.tif" stream = rxr.open_rasterio(outStream) stream = stream[0] stream = stream.where(stream.values!=-32768) # Ignore nodata fig = plt.p(figsize=(20, 10)) ax1 = fig.add_subplot(1, 2, 1) ax2 = fig.add_subplot(1, 2, 2) # 绘制河网图 stream.plot(ax=ax1, cmap='binary', add_colorbar=False) # 定义局部放大的窗口大小、位置（右下角） x_max = stream.coords['x'].max().item() y_min = stream.coords['y'].min().item() res = 30 window_size = 1000 # 在全景图中绘制局部放大区域的范围 rect = mpatches.Rectangle( (x_max-res*window_size, y_min), res*window_size, res*window_size, linewidth=1, edgecolor='b', facecolor='none' ) ax1.add_patch(rect) # 自定义图例：黑色区域为河网，白色区域非河网 patches = [ mpatches.Patch(color='black', label='Is stream'), mpatches.Patch(color='white', label='Is not stream') ] ax1.legend(handles=patches, title='Legend') ax1.set_title('Stream') ax1.set_xticks([]) ax1.set_yticks([]) # 绘制局部放大的河网图 plt.subplot(1,2,2) local_image = stream[-1000:,-1000:] local_image.plot(ax=ax2, cmap='binary', add_colorbar=False) ax2.legend(handles=patches, title='Legend') ax2.set_title('Local Stream') ax2.set_xticks([]) ax2.set_yticks([]) plt.show() del stream del local_image plt.clf()

outStream = temp_dir + "Stream.tif" outSrmLine = temp_dir + "StreamLine.shp" wbt.raster_to_vector_lines( outStream, outSrmLine, callback=my_callback ) outHS = temp_dir + "Hillshade.tif" outSrmLine = temp_dir + "StreamLine.shp" fig, ax = plt.subplots(figsize=(12,12)) hs = rxr.open_rasterio(outHS) hs = hs[0] hs = hs.where(hs.values!=-9999) # Ignore nodata # 将山地阴影图作为底图 streamL = gpd.read_file(outSrmLine) hs.plot(ax=ax, cmap='gray', add_colorbar=False) # 绘制矢量化河网 streamL.plot(ax=ax, legend=True) ax.set_xticks([]) ax.set_yticks([]) ax.set_title('Streams') ax.legend(labels=['Streams'], title='Legend') plt.show() del hs del streamL plt.clf()

outFlowD = temp_dir + "FlowDir.tif" outStream = temp_dir + "Stream.tif" outSrmLink = temp_dir + "StreamLink.tif" wbt.stream_link_identifier( outFlowD, outStream, outSrmLink, esri_pntr=True, zero_background=False, callback=my_callback ) outSrmLink = temp_dir + "StreamLink.tif" link_image = rxr.open_rasterio(outSrmLink) link_image = link_image[0] link_image = link_image.where(link_image.values!=-32768) # Ignore nodata fig = plt.p(figsize=(20, 10)) ax1 = fig.add_subplot(1, 2, 1) ax2 = fig.add_subplot(1, 2, 2) # 绘制河流链图，为不同id的河流链随机分配颜色 link_image.plot(ax=ax1, cmap='Paired', cbar_kwargs={'label': 'ID'}) # 定义局部放大的窗口大小、位置（右下角） x_max = link_image.coords['x'].max().item() y_min = link_image.coords['y'].min().item() res = 30 window_size = 1000 # 在全景图中绘制局部放大区域的范围 rect = mpatches.Rectangle( (x_max-res*window_size, y_min), res*window_size, res*window_size, linewidth=1, edgecolor='b', facecolor='none' ) ax1.add_patch(rect) ax1.set_title('Stream Link') ax1.set_xticks([]) ax1.set_yticks([]) plt.subplot(1,2,2) # 绘制局部放大的河流链图 local_link_image = link_image[-window_size:,-window_size:] local_link_image.plot(ax=ax2, cmap='Paired', cbar_kwargs={'label': 'ID'}) ax2.set_title('Local Stream Link') ax2.set_xticks([]) ax2.set_yticks([]) plt.show() del link_image del local_link_image plt.clf()

outFlowD = temp_dir + "FlowDir.tif" outSrmLink = temp_dir + "StreamLink.tif" outWatershed = output_dir + "WaterShed.tif" wbt.watershed( outFlowD, outSrmLink, outWatershed, esri_pntr=True, callback=my_callback ) outWatershed = output_dir + "WaterShed.tif" outSrmLine = temp_dir + "StreamLine.shp" ws_image = rxr.open_rasterio(outWatershed) ws_image = ws_image[0] ws_image = ws_image.where(ws_image.values!=-32768) # Ignore nodata streamL = gpd.read_file(outSrmLine) # 自定义色带。绘制对象是不同ID的watershed，因此色带颜色应具备以下特征：不重复、各颜色间差异较大 num_colors = len(np.unique(ws_image.values)) spectral_cmap = plt.cm.get_cmap('YlGnBu', num_colors+2) fig, ax = plt.subplots(figsize=(10, 10)) ws_image.plot(cmap=spectral_cmap,cbar_kwargs={'label': 'ID'}) streamL.plot(ax=ax) ax.set_xticks([]) ax.set_yticks([]) ax.set_title('Watershed') plt.show() del ws_image del streamL plt.clf()

outDEM = output_dir + "DEM_P.tif" outMean = temp_dir + "DEM_Mean.tif" windowsize = 200 wbt.mean_filter( outDEM, outMean, filterx=windowsize, filtery=windowsize, callback=my_callback )

import string from itertools import cycle # 直接调用gdal的gdal_calc工具 # expression中的参数为A、B、C...... # input需要按照在expression中与字母一一对应的顺序排列 def rasterCalculator(expression, output, input=[]): os.environ['PROJ_LIB'] = "/opt/conda/share/proj" abs_path = "/opt/conda/lib/python3.9/site-packages/osgeo_utils/gdal_calc.py" if(len(input)==0): return 0 in_str = "" letters_iter = cycle(string.ascii_uppercase) for i in range(len(input)): letter = next(letters_iter) in_str = in_str + " -" + letter +" " + input[i] command = "python " + abs_path + in_str + " --cal=\'" + expression + "\' --outfile=" + output + " --overwrite --NoDataValue=-9999" print(command) # os.system(command) subprocess.run(command, shell=True, stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT) return("success") outDEM = output_dir + "DEM_P.tif" outMean = temp_dir + "DEM_Mean.tif" expression = "A-B" outDiff = temp_dir + "DEM_Diff.tif" rasterCalculator( expression, outDiff, [outDEM,outMean] ) outDiff = temp_dir + "DEM_Diff.tif" diff_image = rxr.open_rasterio(outDiff) diff_image = diff_image[0] diff_image = diff_image.where(diff_image.values!=-9999) # Ignore nodata fig, ax = plt.subplots(figsize=(10, 10)) diff_image.plot(cmap='gray',cbar_kwargs={'label': 'DEM difference'}) ax.set_xticks([]) ax.set_yticks([]) ax.set_title("DEM Difference") plt.show() del diff_image plt.clf()

outDiff = temp_dir + "DEM_Diff.tif" outNeg = output_dir + "Negative.tif" expression = "A<=0" rasterCalculator( expression, outNeg, [outDiff] ) outNeg = output_dir + "Negative.tif" neg_image = rxr.open_rasterio(outNeg) neg_image = neg_image[0] neg_image = neg_image.where(neg_image.values!=-9999) # Ignore nodata fig, ax = plt.subplots(figsize=(10, 10)) neg_image.plot(cmap='binary', add_colorbar=False) patches = [ mpatches.Patch(color='black', label='Is negative'), mpatches.Patch(color='white', label='Is not negative') ] ax.legend(handles=patches, title='Legend') ax.set_xticks([]) ax.set_yticks([]) ax.set_title("Negative") plt.show() del neg_image plt.clf()

# 调用大平台模型库中的矢量转栅格工具，修改了部分代码：指定field def pydde_Vector2Raster(shapefile_path, output_raster_path, pixel_size, field = "class", extent: tuple = None): input_shp = ogr.Open(shapefile_path) # getting layer information of shapefile. shp_layer = input_shp.GetLayer() # get extent values to set size of output raster. if extent is None: x_min, x_max, y_min, y_max = shp_layer.GetExtent() else: x_min, x_max, y_min, y_max = extent # calculate size/resolution of the raster. x_res = int((x_max - x_min) / pixel_size) y_res = int((y_max - y_min) / pixel_size) # print(x_res, y_res) # get GeoTiff driver by image_type = 'GTiff' driver = gdal.GetDriverByName(image_type) # passing the filename, x and y direction resolution, no. of bands, new raster. new_raster = driver.Create(output_raster_path, x_res, y_res, 1, gdal.GDT_Byte) # transforms between pixel raster space to projection coordinate space. new_raster.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size)) # get required raster band. band = new_raster.GetRasterBand(1) # assign no data value to empty cells. no_data_value = -99999 band.SetNoDataValue(no_data_value) band.FlushCache() # Core conversion method gdal.RasterizeLayer(new_raster, [1], shp_layer, burn_values=[255],options=[f"ATTRIBUTE={field}"]) # adding a spatial reference spatialRef = shp_layer.GetSpatialRef() # print(spatialRef) if spatialRef is None: new_rasterSRS = osr.SpatialReference() new_rasterSRS.ImportFromEPSG(4326) new_raster.SetProjection(new_rasterSRS.ExportToWkt()) else: new_raster.SetProjection(spatialRef.ExportToWkt()) print("Conversion successful") def getRasterExtent(raster_path): ds = gdal.Open(raster_path) width = ds.RasterXSize height = ds.RasterYSize trans = ds.GetGeoTransform() x_min = trans[0] y_max = trans[3] x_max = trans[0] + width * trans[1] + height * trans[2] y_min = trans[3] - abs(width * trans[4]) - abs(height * trans[5]) res = trans[1] return (x_min, x_max, y_min, y_max) def getRasterResolution(raster_path): ds = gdal.Open(raster_path) trans = ds.GetGeoTransform() return trans[1] outDEM = output_dir + "DEM_P.tif" outGBLU = output_dir + "GBLU_P.shp" outVec2ras = temp_dir + "GBLU_ras_P.tif" field = "gridcode" ext = getRasterExtent(outDEM) #面转栅格的范围与DEM保持一致 pixel_size = getRasterResolution(outDEM) pydde_Vector2Raster(outGBLU, outVec2ras, pixel_size, field, extent=ext) outVec2ras = temp_dir + "GBLU_ras_P.tif" class_image = rxr.open_rasterio(outVec2ras) class_image = class_image[0] class_image = class_image.where(class_image.values!=0) # Ignore nodata landform_type = { "Low altitude plain":101, "Middle altitude plain":102, "Low altitude hill":211, "Middle altitude hill":212, "Middle altitude low-relief mountain":222 } # 每个code对应的地貌类型 fig, ax = plt.subplots(figsize=(10, 10)) # 自定义色带，按照code对地貌类型设色 num_colors = len(landform_type)+1 cmap = plt.cm.get_cmap('YlGn_r', num_colors) landform_cmap = ListedColormap([cmap(i) for i in range(num_colors)]) boundaries = list(landform_type.values()) norm = BoundaryNorm(boundaries, cmap.N, extend='max') class_image.plot(cmap=landform_cmap, levels=boundaries+[230], add_colorbar=False) # 自定义图例 for value in boundaries: plt.scatter([], [], color=landform_cmap(norm(value)-1))#, label=str(value) ax.legend(title="Landform Type", loc="upper right", labels=list(landform_type.keys())) ax.set_title("Basic Landform Unit") ax.set_xticks([]) ax.set_yticks([]) plt.show() del class_image

outVec2ras = temp_dir + "GBLU_ras_P.tif" outPlain = temp_dir + "Plain_ras.tif" reclass_vals='1;0;200;0;200;1000' #重分类表达式：0-200->1 ; 200-1000->0 wbt.reclass( outVec2ras, outPlain, reclass_vals, assign_mode=False, callback=my_callback ) outPlain = temp_dir + "Plain_ras.tif" plain_image = rxr.open_rasterio(outPlain) plain_image = plain_image[0] plain_image = plain_image.where(plain_image.values!=-32768) # Ignore nodata fig, ax = plt.subplots(figsize=(10, 10)) plain_image.plot(cmap='binary', add_colorbar=False) # 自定义图例：黑色为平原区域，白色非平原区域 patches = [ mpatches.Patch(color='black', label='Is plain'), mpatches.Patch(color='white', label='Is not plain') ] ax.legend(handles=patches, title='Legend') ax.set_xticks([]) ax.set_yticks([]) ax.set_title("Plain Area") plt.show() del plain_image plt.clf()

outNeg = output_dir + "Negative.tif" outPlain = temp_dir + "Plain_ras.tif" outBottom = output_dir + "Bottom.tif" expression = "A*B" rasterCalculator( expression, outBottom, [outNeg,outPlain] ) outBottom = output_dir + "Bottom.tif" bottom_image = rxr.open_rasterio(outBottom) bottom_image = bottom_image[0] bottom_image = bottom_image.where(bottom_image.values!=-9999) # Ignore nodata fig, ax = plt.subplots(figsize=(10, 10)) bottom_image.plot(cmap='binary', add_colorbar=False) # 自定义图例：黑色为沟底区域，白色非沟底区域 patches = [ mpatches.Patch(color='black', label='Is bottom'), mpatches.Patch(color='white', label='Is not bottom') ] ax.legend(handles=patches, title='Legend') ax.set_xticks([]) ax.set_yticks([]) ax.set_title("Bottom of Trench") plt.show() del bottom_image plt.clf()

outNeg = output_dir + "Negative.tif" outBottom = output_dir + "Bottom.tif" outSlope = output_dir +"Slope.tif" expression = "A-B" rasterCalculator( expression, outSlope, [outNeg,outBottom] ) outSlope = output_dir +"Slope.tif" slope_image = rxr.open_rasterio(outSlope) slope_image = slope_image[0] slope_image = slope_image.where(slope_image.values!=-9999) # Ignore nodata fig, ax = plt.subplots(figsize=(10, 10)) slope_image.plot(cmap='binary', add_colorbar=False) # 自定义图例：黑色为沟坡区域，白色非沟坡区域 patches = [ mpatches.Patch(color='black', label='Is slope'), mpatches.Patch(color='white', label='Is not slope') ] ax.legend(handles=patches, title='Legend') ax.set_xticks([]) ax.set_yticks([]) ax.set_title("Slope of Trench") plt.show() del slope_image plt.clf()

outSlope = output_dir +"Slope.tif" features = outWatershed = output_dir + "WaterShed.tif" outSlopeArea = temp_dir + "SlopeArea.tif" wbt.zonal_statistics( outSlope, features, outSlopeArea, stat="total", out_table=None, callback=my_callback )

outWatershed = output_dir + "WaterShed.tif" outWsArea = temp_dir + "WaterShedArea.tif" wbt.raster_area( outWatershed, outWsArea, out_text=False, units="grid cells", zero_back=False, callback=my_callback )

outSlopeArea = temp_dir + "SlopeArea.tif" outWsArea = temp_dir + "WaterShedArea.tif" outStat = output_dir + "SlopeRatio.tif" expression = "A/B" rasterCalculator( expression, outStat, [outSlopeArea,outWsArea] ) outStat = output_dir + "SlopeRatio.tif" outSrmLine = temp_dir + "StreamLine.shp" result = rxr.open_rasterio(outStat) result = result[0] result = result.where(result.values!=-9999) # Ignore nodata streamL = gpd.read_file(outSrmLine) fig, ax = plt.subplots(figsize=(12,12)) cmap = plt.cm.get_cmap('RdYlBu_r') # boundaries = np.linspace(0, 1, num=6) # norm = BoundaryNorm(boundaries, cmap.N) result.plot(cmap=cmap, add_colorbar=True, cbar_kwargs={'label': 'ratio'}) streamL.plot(ax=ax) ax.set_xticks([]) ax.set_yticks([]) ax.set_title("Slope Ratio") plt.show()