博客

1.利用xarray导入tif格式的DEM影像,并让其可视化。

参考博客：

主要文章

和

博客1

及

博客2

。

代码如下：

首先导入相关的包：

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from osgeo import gdal
import cmaps
import matplotlib as mpl
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib import rcParams

其次导入文件：

BJ_tiff=r'D:\google\earth_engine\Gee_Gis_file\GuiZhou\毕节市\Bijie_Tif\Bijie.tif'
BJ_DEM=r'D:\google\earth_engine\Gee_Gis_file\DEM\BJ_DEM_Elevation.tif'
HeNan_tiff=r'D:\google\earth_engine\Gee_Gis_file\DEM\HeNan_Elevation.tif'

然后读取文件并可视化显示：

if __name__=='__main__':
   BJ_dem=xr.open_rasterio(BJ_DEM)
   X, Y = np.meshgrid(BJ_dem.x,BJ_dem.y)
   Z=np.array(BJ_dem.sel(band=1))
   proj=ccrs.PlateCarree()
   extent = [min(BJ_dem.x),max(BJ_dem.x),min(BJ_dem.y),max(BJ_dem.y)]
   fig = plt.figure(figsize=(14, 10),dpi=600)
   ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
   ax.set_xticks(np.arange(extent[0], extent[1] + 1,0.3), crs = proj)
   ax.set_yticks(np.arange(extent[-2], extent[-1] + 1,0.3), crs = proj)
   ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
   ax.yaxis.set_major_formatter(LatitudeFormatter())
   ax.set_extent(extent, crs=ccrs.PlateCarree())#增加经度和纬度
   shp_path = r'D:\google\earth_engine\Gee_Gis_file\GuiZhou\Bijie_IL\Bijie_IL.shp'
   proj = ccrs.PlateCarree() 
   reader = Reader(shp_path)
   provinces = cfeature.ShapelyFeature(reader.geometries(), proj,edgecolor='k', facecolor='none',alpha=1)
   ax.add_feature(provinces, linewidth=0.65)
   lev=np.arange(0,2400,200)
   ticks=[0,200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400]
   cf=ax.contourf(X,Y,Z,levels=lev,extend='both',transform=ccrs.PlateCarree(),cmap=cmaps.MPL_terrain)
   b=plt.colorbar(cf,ticks=ticks,shrink=0.55,orientation='vertical',extend='both',pad=0.015,aspect=20)
# plt.savefig('F:/Rpython/lp32/plot11.3.png',dpi=600)
   plt.show()

简单可视化为：