python绘制中国版dem地图

背景

前几天看到一个人分享的中国的dem数据可视化实在是太好看了。这里我想用python复刻一下。

过程

1. 数据：

数据已经上传到我的网盘上了，关注公众号【world of statistics】然后回复【气象数据】即可获得。

2. 代码

# 导入包
import rasterio
import numpy as np
import matplotlib.pyplot as plt
from getchinamap.getchinamap import DownloadChmap
import rasterio.mask
import cartopy.crs as ccrs
import matplotlib as mpl
import datetime
dl_engine = DownloadChmap()

# 获得数据
dem_raw_data = rasterio.open("数据集/ChinaDEM1km_WGS84.tif")

这个数据的详情如下：

# 这里以西藏自治区为例：

# china_boundary_gdf = dl_engine.download_country(target='边界')
province_name = '西藏自治区'
china_boundary_gdf = dl_engine.download_province(province_name=province_name,target='边界')
china_boundary_gdf_sub = dl_engine.download_province(province_name=province_name,target='市')
# china_boundary_gdf = dl_engine.download_country(target='边界')
# china_boundary_gdf_sub = dl_engine.download_country(target='省')
china_boundary_gdf_sub

# 对dem做裁切
out_image, out_transform = rasterio.mask.mask(dem_raw_data, china_boundary_gdf['geometry'], crop=True)
out_meta = dem_raw_data.meta
out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})
                 
                 
clean_data = out_image.copy().astype('float')
clean_data = clean_data[0, :, :]
clean_data[clean_data == out_meta['nodata']] = np.nan

3. 使用`ccrs.PlateCarree()`投影

接下来是使用ccrs.PlateCarree()投影，使用这个投影的原因是可视化快。代码如下：


mat_extent = [getattr(out_transform, 'c'),
              getattr(out_transform, 'c') + clean_data.shape[1] * getattr(out_transform, 'a'),
              getattr(out_transform, 'f') + clean_data.shape[0] * getattr(out_transform, 'e'),
              getattr(out_transform, 'f')]

colors = ["#33A02C", "#B2DF8A", "#FDBF6F", "#1F78B4", "#999999", "#E31A1C", "#E6E6E6"


    
, "#A6CEE3"]

myccrs = ccrs.PlateCarree()

fig, ax = plt.subplots(figsize=(10, 5), dpi=150, subplot_kw={'projection': myccrs})

ax.gridlines(draw_labels=True,
             linewidth=2, color='gray', alpha=0.5, linestyle='--')

ax_im_bar = ax.imshow(clean_data, origin='upper', extent=mat_extent, transform=myccrs,
                      cmap=mpl.colors.LinearSegmentedColormap.from_list("mypalette", colors, N=1000))

china_boundary_gdf_sub.boundary.plot(ax=ax, color='black', linewidth=1)

fig.colorbar(ax_im_bar, orientation='vertical')
# plt.tight_layout()
#
# ax.set_title("china dem")
# ax.set_ylabel("y label")
# ax.set_xlabel("x label")

4. 使用`Orthographic`投影

但是上面的图不太好看，这里使用Orthographic投影再画一次。但是这个渲染的速度非常慢，需要等待几分钟⌛️～

longitude = np.array(
    [getattr(out_transform, 'c') + i * getattr(out_transform, 'a') for i in range(clean_data.shape[1])])
latitude = np.array(
    [getattr(out_transform, 'f') + i * getattr(out_transform, 'e') for i in range(clean_data.shape[0])])

colors = ["#33A02C", "#B2DF8A", "#FDBF6F", "#1F78B4", "#999999", "#E31A1C", "#E6E6E6", "#A6CEE3"]
center_china_point = china_boundary_gdf.centroid[0]

# 新的兰伯特投影
new_crs = ccrs.Orthographic(central_longitude=center_china_point.x,
                                central_latitude=center_china_point.y)
myccrs = ccrs.PlateCarree()

fig, ax = plt.subplots(figsize=(10, 10), dpi=100, subplot_kw={'projection': new_crs})

ax.gridlines(draw_labels=True,
             linewidth=2, color='gray', alpha=0.5, linestyle='--')
ax_im_bar = ax.contourf(longitude, latitude, clean_data,  # [::-1,:]
                            np.arange(0, np.nanmax(clean_data), 100.0),
                            transform=ccrs.PlateCarree(),
                            cmap=mpl.colors.LinearSegmentedColormap.from_list("mypalette", colors, N=1000),
                            extend='both', origin='upper')
china_boundary_gdf_sub.to_crs(new_crs.proj4_init).boundary.plot(ax=ax, edgecolor='black', linewidth=1)
fig.colorbar(ax_im_bar, orientation='vertical')
# fig.savefig('test1226.png')
fig.savefig(f"test1226_china_{datetime.datetime.now().strftime('%Y-%m-%d %X')}.png")