进阶！dask解决超高精度tif读取与绘图难问题

用户11172986

发布于 2024-06-20 18:05:21

740

发布于 2024-06-20 18:05:21

文章被收录于专栏：气python风雨气python风雨

前言

又有读者来信要求如下：希望小编帮忙看看能不能解决。是关于能不能在已经截取出来的省份中添加对应的dem地形呢，并且根据需要添加上需要的城市所在的地理位置，比如在已绘制的图中标注出三亚的所在地

数据：地形tif文件难点：文件格点过多，可视化会爆内存解决办法：dask延迟加载，分块读取，绘图方式采用imshow 镜像：气象分析3.9

In [1]:

!pip install rioxarray -i https://pypi.mirrors.ustc.edu.cn/simple/
!pip install cnmaps -i https://pypi.mirrors.ustc.edu.cn/simple/

原程序

In [6]:

"""
Created on Tue Jan 23 20:13:17 2024

@author: Ritchie
"""

import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch 
from scipy.interpolate import Rbf
import shapefile
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib import colors 
import matplotlib as mpl
from matplotlib.axes import Axes

proj = ccrs.PlateCarree()  # 地图投影：圆柱投影
fig, ax = plt.subplots(figsize = (15, 15), subplot_kw=dict(projection=proj))  # 画布
#----------------------------------------------------------------------------

shpname = '/home/mw/input/dem5930/求助小编'
shpname1 = 'F:/Study/中国基础数据/'
#fig, ax = plt.subplots(figsize = (15, 15), subplot_kw=dict(projection=proj))  # 画布
ax.set_extent([108.3, 111, 18, 20.2], ccrs.PlateCarree())
#ax.gridlines(linestyle='-')
ax.add_geometries(shpreader.Reader(os.path.join(shpname, 'Hainan_province.shp')).geometries(),ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1)#地理信息
ax.add_geometries(shpreader.Reader(os.path.join(shpname, 'Hainan_city.shp')).geometries(),ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1)
#ax.add_geometries(shpreader.Reader(os.path.join(shpname, 'Hainan_county.shp')).geometries(),ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1)
#----------------------------------------------------------------------------------------------

#----------------------------------------------------------------------------
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())

ax.set_xticks(np.arange(108.3, 111.6, .3),crs = ccrs.PlateCarree())    # x轴
ax.set_yticks(np.arange(18, 20.3, .2),crs = ccrs.PlateCarree())      # y轴
ax.tick_params(labelsize=17)

shp绘制在气象分析3.7环境下是正常绘制，有城市划分，但使用气象分析3.9会出现以上不全的情况，笔者无法解决

实现地形倒也不难，之前画过很多关于地形的图例如Python 绘制山体阴影+雷达图

导入库

In [2]:

import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch 
from scipy.interpolate import Rbf
import shapefile
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib import colors 
import matplotlib as mpl
from matplotlib.axes import Axes
import rioxarray
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import dask.array as da
from cnmaps import get_adm_maps, draw_maps, clip_contours_by_map

数据读取

In [3]:

# 读取地形tif文件
data = rioxarray.open_rasterio("/home/mw/input/dem5930/海南省WGS84.tif")
data

地形数据读取成功，但是在实际绘图时常常会爆内存，怎么回事一看地形数据是481805534 values with dtype=int16 那没事了这时候就需要dask出动

什么是dask

Dask 是一个灵活的并行计算库，旨在处理大型数据集。它提供了一种能够处理比内存更大的数据集的方法，并能够以并行和延迟加载的方式执行计算任务。

主要特点包括：

并行化: Dask 可以自动并行执行多个任务，从而充分利用多核 CPU 或者集群资源来加速计算。

延迟加载: Dask 支持延迟加载（lazy evaluation），这意味着它只有在真正需要执行计算时才会加载数据并执行操作。

分布式计算: Dask 支持分布式计算，可以在分布式环境中运行，处理跨多台计算机的大规模数据集。

适用范围: Dask 可以用于各种数据类型，包括数组、DataFrame 和机器学习模型等。

总之，Dask 提供了一种便捷的方式来处理大型数据集，并且能够有效地进行并行计算，从而加速数据处理过程。

In [4]:

# 读取地形tif文件（使用延迟加载）
data1 = rioxarray.open_rasterio("/home/mw/input/dem5930/海南省WGS84.tif", chunks={'band': 1, 'x': 1024, 'y': 1024})
data1

绘图部分

In [10]:

import cmaps

cmap = cmaps.MPL_terrain

In [6]:

proj = ccrs.PlateCarree()  # 地图投影：圆柱投影
fig, ax = plt.subplots(figsize = (15, 10), subplot_kw=dict(projection=proj))  # 画布
#----------------------------------------------------------------------------

shpname = '/home/mw/input/dem5930/求助小编'
ax.set_extent([data1.x.min(), data1.x.max(), data1.y.min(), data1.y.max()], ccrs.PlateCarree())

ax.add_geometries(shpreader.Reader(os.path.join(shpname, 'Hainan_province.shp')).geometries(),ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1)#地理信息
ax.add_geometries(shpreader.Reader(os.path.join(shpname, 'Hainan_city.shp')).geometries(),ccrs.PlateCarree(), facecolor='none', edgecolor='w', linewidth=1)
ax.add_geometries(shpreader.Reader(os.path.join(shpname, 'Hainan_county.shp')).geometries(),ccrs.PlateCarree(), facecolor='none', edgecolor='w', linewidth=1)
#----------------------------------------------------------------------------------------------

#----------------------------------------------------------------------------
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())

ax.set_xticks(np.arange(108.6, 111.6, .3),crs = ccrs.PlateCarree())    # x轴
ax.set_yticks(np.arange(18, 20.1, .2),crs = ccrs.PlateCarree())      # y轴
ax.tick_params(labelsize=17)


# 使用 Dask array 进行绘图
chunked_data = da.from_array(data1.values, chunks=(1, data1.shape[1]//10, data1.shape[2]//10))
# 绘制地形图
img =ax.imshow(chunked_data[0], extent=(data1.x.min(), data1.x.max(), data1.y.min(), data1.y.max()), cmap=cmap)

ax.add_geometries(shpreader.Reader(os.path.join(shpname, 'Hainan_province.shp')).geometries(),ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1)#地理信息
ax.add_geometries(shpreader.Reader(os.path.join(shpname, 'Hainan_city.shp')).geometries(),ccrs.PlateCarree(), facecolor='none', edgecolor='w', linewidth=1)
ax.add_geometries(shpreader.Reader(os.path.join(shpname, 'Hainan_county.shp')).geometries(),ccrs.PlateCarree(), facecolor='none', edgecolor='w', linewidth=1)

# 标记三亚的中心
ax.text(109.518, 18.257, '三亚',color='w', fontsize=18, ha='center', transform=ccrs.PlateCarree())

# 添加色卡
cbar = fig.colorbar(img, ax=ax, orientation='vertical')
cbar.set_label('Elevation')

plt.show()

cnmaps优化版（不读取shp）

In [6]:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from cnmaps import get_adm_maps, draw_maps

fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())

# 使用 Dask array 进行绘图
chunked_data = da.from_array(data1.values, chunks=(1, data1.shape[1]//10, data1.shape[2]//10))
# 绘制地形图
img =ax.imshow(chunked_data[0], extent=(data1.x.min(), data1.x.max(), data1.y.min(), data1.y.max()), cmap=cmap)


# 标记三亚的中心
ax.text(109.518, 18.257, '三亚',color='w', fontsize=18, ha='center', transform=ccrs.PlateCarree())

# 添加色卡
cbar = fig.colorbar(img, ax=ax, orientation='vertical')
cbar.set_label('Elevation')

province = get_adm_maps(province='海南省', record='first', only_polygon=True)
draw_maps(get_adm_maps(province='海南省', level='市'), color='w', linewidth=0.8)
ax.set_extent([data1.x.min(), data1.x.max(), data1.y.min(), data1.y.max()], ccrs.PlateCarree())
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())

ax.set_xticks(np.arange(108.6, 111.6, .3),crs = ccrs.PlateCarree())    # x轴
ax.set_yticks(np.arange(18, 20.1, .2),crs = ccrs.PlateCarree())      # y轴
ax.tick_params(labelsize=17)

plt.show()