[python] 画图：cartopy笔记（一）基础知识与误差分布图

已于 2023-06-27 19:28:08 修改
阅读量1k 收藏 7
点赞数 1
分类专栏： 笔记 文章标签： python
于 2022-04-06 18:12:34 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_52430372/article/details/123984810
版权

参考链接：

基础知识理解：

摸鱼的气象&Python：绘图基础和地图投影_哔哩哔哩_bilibili

对应代码：摸鱼的气象&Python - Heywhale.com

《气py：Python气象资料处理与可视化》李开元、豆京华

​​​​​​Cartopy 系列：从入门到放弃 - 炸鸡人博客 (zhajiman.github.io)

画图实例：

​​​​​​(29条消息) 小白学习cartopy画地图的第一天（中国行政区域图，含南海）_野生的气象小流星的博客-CSDN博客_cartopy绘制中国省份图

​​​​​​(29条消息) 小白学习cartopy气象画地图的第二天（中国区域，陆地温度分布图）_野生的气象小流星的博客-CSDN博客_cartopy 气象

最终效果图如下：

目录

参考链接：

1.cartopy简介

2.画图流程及常用函数

2.1.创建画布、绘画区（figure、axes）

Figure：画布

Axes：绘图区

2.2 画地图

投影方式

3.画图实例：华南地区降水量mae分布图

1.cartopy简介

cartopy是一个为地理空间数据处理而设计的python包，用于生成地图和其他地理空间数据分析。本身不具有可视化功能，但可以为Matplotlib提供地理坐标系的转换，实现气象专业图形的绘制。

2.画图流程及常用函数

2.1.创建画布、绘画区（figure、axes）

import matplotlib.pyplot as plt

Figure：画布

Figure(figsize=None, dpi=None, facecolor=None, edgecolor=None, linewidth=0.0)

图像尺寸、像素点、背景颜色、边框颜色、边框宽度

import matplotlib.pyplot as plt

fig=plt.figure(figsize=(4,3),dpi=200,facecolor='blue',edgecolor='black',linewidth=2)

plt.show()

Axes：绘图区

figure对象有两种方法建立axes: plt.figure().add_axes（）、 plt.figure().add_subplot（）

其中，add_subplots属于axes的特殊情况

1).add_axes

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_axes((0,0,3,2))#以(0,0)点为左下的起点坐标，选择一个长为3，宽为2的绘图区 
plt.show()

2).add_subplot

import matplotlib.pyplot as plt
fig = plt.figure()
ax = add_subplot(2，2，1)
plt.show()

3).GeoAxes：地图投影，把坐标投影为地理坐标

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig = plt.figure()
ax = fig.add_axes([0,0,3,2],projection = ccrs.PlateCarree())
plt.show()

2.2 画地图

投影方式

有很多，最常用等距圆柱投影，对数据没有变形

可参考：摸鱼的气象&Python：绘图基础和地图投影_哔哩哔哩_bilibili

1).等距圆柱投影   ccrs.PlateCarree(central_longitude=0.0)

2).正形兰伯特投影   ccrs.LambertConformal( central_longitude=-96.0, central_latitude=39.0,
    standard_parallels=None, cutoff=-30,)

中心经度、中心维度、区间（内部数据为等比例投影）、截止纬度

填色（contour、contourf）

matplotlib.axes.Axes.contour:绘制等值线

matplotlib.axes.Axes.contourf：给等值线填色

ax.contourf(lon.lat,z,levels,cmap,transform,zorder)

levels:整数n或数组，n:画n条线，数组：按数组值画等值线

cmap:当color未被指定时，使用自己指定的cmap或官方默认的colormap方案。

fig = plt.figure(figsize=(12,8))
proj = ccrs.PlateCarree(central_longitude=180)
ax = fig.add_axes([0.1, 0.1, 0.8, 0.6],projection = proj)
c1 = ax.contour(lon,lat,z,levels=np.arange(5400,6000,50),
 cmap=plt.cm.bwr,transform=ccrs.PlateCarree())

c2 = ax.contourf(lon,lat,z,levels=np.arange(5400,6000,50),
extend='max' ,cmap=plt.cm.bwr,
transform=ccrs.PlateCarree(),zorder=0)

增加（刻度线、标题）等

可参考：摸鱼的气象&Python：几个必不可少的地理绘图函数_哔哩哔哩_bilibili

import cartopy.crs as ccrs
import cartopy.mpl.ticker as cticker
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(12,8))
proj = ccrs.PlateCarree(central_longitude=180)
ax = fig.add_axes([0.1, 0.1, 0.8, 0.6],projection = proj)

ax.set_xticks(np.arange(leftlon,rightlon+60,60), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(lowerlat,upperlat+30,30), crs=ccrs.PlateCarree())
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

3.画图实例：华南地区降水量mae分布图

导入库

#参考链接：
# https://blog.csdn.net/weixin_42372313/article/details/113665724?spm=1001.2014.3001.5502
# https://zhajiman.github.io/post/cartopy_introduction/#%E7%AE%80%E4%BB%8B

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter

一、加载数据mae、中国地图文件

plt.rcParams["font.sans-serif"] = ["DengXian"]  # 用来正常显示中文标签
plt.rc('axes', unicode_minus=False)  # 解决负号显示问题

# 加载mae mae（step,latitude,longitude）,取超前0天对应的mae
correlation = mae.sel(step = '0 days')
lon = mae.longitude
lat = mae.latitude

#画地图底图(网址是：*https://gmt-china.org/data/*下载里面的CN-border-La.dat)
## 这一步是读取CN - border - La.da文件，目的是读取里面每个点的经纬度，用于给底图加上行政边界

 with open('CN-border-La.dat') as src:
        context = src.read()
        blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
        borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]

二、画图

#1、创建画布
fig = plt.figure(figsize=[8, 8])#画布
ax = plt.axes(projection=ccrs.PlateCarree())#ax:绘图区，projection:投影方式为PlateCarree


#2、画底图

# 这里就用到了当时上面文件里面的坐标点了，其实就是用点画线，
# line[0::2]就是X，line[1::2]就是Y，就是用xy画平面图，一组是经度，一组是纬度。
    
 for line in borders:
        ax.plot(line[0::2], line[1::2], '-', lw=1.5, color='k',transform=ccrs.Geodetic())


# 咱们画出来的是世界地图，这个来框处要显示的区域
xmin, xmax = lon.min(), lon.max()
ymin, ymax = lat.min(), lat.max()
ax.set_extent([xmin, xmax, ymin,ymax],crs=ccrs.PlateCarree())


#3、画图填色
#使用contourf进行填色，下面先设置它的参数

#设置色标
import math
t_max = math.ceil(correlation.max().values)
t_min = math.floor(correlation.min().values)
n = (t_max - t_min)/10
cbar_kwargs = {
    'orientation': 'vertical',
    'label': '降水量误差(mm)',
    'shrink': 0.8,
    'ticks': np.arange(t_min, t_max, n),
    'pad': 0.05,
    'shrink': 0.65
}

# 设置画数据的精度，这里设置了从数据的最小值到最大值，按0.1mm画降水量
levels = np.arange(t_min, t_max, 0.1)

# 本质上是一个画等值线并填色的函数，被cartopy封装了一下，参数的含义可以去查#matplotlib.pyplot.contourf
# 参数含义：levels设置等值线密度，cmap颜色，cbar_kwargs色标，transform变换
correlation.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r',
                        cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())


#4、画经纬度（即网格线）、刻度线、标题
#设置x、y刻度
#以5°为间隔计算x、y坐标刻度值，使用PlateCarree投影方式换算坐标，得到x、y的地理坐标
ax.set_xticks(
    np.arange(xmin, xmax + 1, np.round((xmax + 1 - xmin) / 5)),
    crs=ccrs.PlateCarree(),
)
ax.set_yticks(
    np.arange(ymin, ymax + 1, np.round((ymax + 1 - ymin) / 5)),
    crs=ccrs.PlateCarree(),
)    

# 设置刻度格式为经纬度格式
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())


#设置标题
ax.set_title('华南地区总降水量mae分布图', color='blue', fontsize=20)

三、显示图形

plt.show()

四、完整代码

#参考链接：
# https://blog.csdn.net/weixin_42372313/article/details/113665724?spm=1001.2014.3001.5502
# https://zhajiman.github.io/post/cartopy_introduction/#%E7%AE%80%E4%BB%8B

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter

class Error:
    def __init__(self,y_pred2,y_obs):
        self.y_pred2 = y_pred2
        self.y_obs = y_obs

    def Rmse(self):
        '''
        计算订正值与真实值之间的rmse
        :param self.y_pred2: 二次预测值（订正）
        :param self.y_obs: 观测值
        :return: rmse(i,j,l)、rmse_mean(l)
        '''
        err = self.y_pred2 - self.y_obs
        err = err**2
        err = err.mean(dim = 'time')
        rmse = err**0.5
        rmse_mean = rmse.mean(dim=['latitude','longitude'])
        return rmse,rmse_mean

    def Mae(self):
        '''
        计算订正值与真实值之间的rmse
        :param self.y_pred2: 二次预测值（订正）
        :param self.y_obs: 观测值
        :return: mae(i,j,l)、mae_mean(l)
        '''
        err = self.y_pred2 - self.y_obs
        err = abs(err)
        mae = err.mean(dim = 'time')
        mae_mean = mae.mean(dim=['latitude','longitude'])
        return mae,mae_mean

    def Tcc(self):
        obs = self.y_obs - self.y_obs.mean("step")
        pred = self.y_pred2 - self.y_pred2.mean("step")

        cov = (obs * pred).sum("step")
        var_var = np.sqrt((obs ** 2).sum("step") * (pred ** 2).sum("step"))
        tcc = cov / var_var
        tcc_mean = tcc.mean('time')
        return tcc,tcc_mean

    def Acc(self):
        self.y_obs = self.y_obs.stack(pos=("latitude", "longitude"))  # 拉成一个维度
        self.y_pred2 = self.y_pred2.stack(pos=("latitude", "longitude"))  # 拉成一个维度
        obs = self.y_obs - self.y_obs.mean("pos")
        pred = self.y_pred2 - self.y_pred2.mean("pos")

        cov = (obs * pred).sum("pos", skipna=True)
        var_var = np.sqrt((obs ** 2).sum("pos", skipna=True) * (pred ** 2).sum("pos", skipna=True))
        acc = cov / var_var
        acc_mean = acc.mean('time')
        return acc,acc_mean

plt.rcParams["font.sans-serif"] = ["DengXian"]  # 用来正常显示中文标签
plt.rc('axes', unicode_minus=False)  # 解决负号显示问题

if __name__ == '__main__':
    # 一、加载数据mae、中国地图文件
    path2 = '/home/gyy/gyy-project/model_set/pred2/source_pred2.nc'  # y_pred2
    path1 = '/home/gyy/gyy-project/model_set/pred2/source_test.nc'  # y_obs
    path0 = '/home/gyy/gyy-project/model_set/pred2/x_mean49.nc'  # y_pred1
    y_pred2 = xr.open_dataarray(path2)
    y_obs = xr.open_dataarray(path1)
    err = Error(y_pred2,y_obs)
    mae,mae_mean = err.Mae()

    # 加载mae mae（step,latitude,longitude）,取超前0天对应的mae
    correlation = mae.sel(step='0 days')
    lon = mae.longitude
    lat = mae.latitude

    # 画地图底图(网址是：*https://gmt-china.org/data/*下载里面的CN-border-La.dat)
    ## 这一步是读取CN - border - La.da文件，目的是读取里面每个点的经纬度，用于给底图加上行政边界

    with open('CN-border-La.dat') as src:
        context = src.read()
        blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
        borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]

    # 二、画图
    # 1、创建画布
    fig = plt.figure(figsize=[8, 8])  # 画布
    ax = plt.axes(projection=ccrs.PlateCarree())  # ax:绘图区，projection:投影方式为PlateCarree

    # 2、画底图

    # 这里就用到了当时上面文件里面的坐标点了，其实就是用点画线，
    # line[0::2]就是X，line[1::2]就是Y，就是用xy画平面图，一组是经度，一组是纬度。

    for line in borders:
        ax.plot(line[0::2], line[1::2], '-', lw=1.5, color='k', transform=ccrs.Geodetic())

    # 咱们画出来的是世界地图，这个来框处要显示的区域
    xmin, xmax = lon.min(), lon.max()
    ymin, ymax = lat.min(), lat.max()
    ax.set_extent([xmin, xmax, ymin, ymax], crs=ccrs.PlateCarree())

    # 3、画图填色
    # 使用contourf进行填色，下面先设置它的参数

    # 设置色标
    import math

    t_max = math.ceil(correlation.max().values)
    t_min = math.floor(correlation.min().values)
    n = (t_max - t_min) / 10
    cbar_kwargs = {
        'orientation': 'vertical',
        'label': '降水量误差(mm)',
        'shrink': 0.8,
        'ticks': np.arange(t_min, t_max, n),
        'pad': 0.05,
        'shrink': 0.65
    }

    # 设置画数据的精度，这里设置了从数据的最小值到最大值，按0.1mm画降水量
    levels = np.arange(t_min, t_max, 0.1)

    # 本质上是一个画等值线并填色的函数，被cartopy封装了一下，参数的含义可以去查#matplotlib.pyplot.contourf
    # 参数含义：levels设置等值线密度，cmap颜色，cbar_kwargs色标，transform变换
    correlation.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r',
                              cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())

    # 4、画经纬度（即网格线）、刻度线、标题
    # 设置x、y刻度
    # 以5°为间隔计算x、y坐标刻度值，使用PlateCarree投影方式换算坐标，得到x、y的地理坐标
    ax.set_xticks(
        np.arange(xmin, xmax + 1, np.round((xmax + 1 - xmin) / 5)),
        crs=ccrs.PlateCarree(),
    )
    ax.set_yticks(
        np.arange(ymin, ymax + 1, np.round((ymax + 1 - ymin) / 5)),
        crs=ccrs.PlateCarree(),
    )

    # 设置刻度格式为经纬度格式
    ax.xaxis.set_major_formatter(LongitudeFormatter())
    ax.yaxis.set_major_formatter(LatitudeFormatter())

    # 设置标题
    ax.set_title('华南地区总降水量mae分布图', color='blue', fontsize=20)

    # 三、显示图形
    plt.show()
白栈央
关注
  • 1
    点赞
  • 踩
  • 7
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
Python绘制地理图--Cartopy基础
qq_40497403的博客
07-30 1万+
常用地图底图的绘制一般由Basemap或者cartopy模块完成，由于Basemap库是基于python2开发的一个模块，目前已经不开发维护。故简单介绍cartopy模块的一些基础操作。 将地球三维球体投影到二维面上，减少失真。主要方式有默认投影（PlateCarree）、兰勃脱投影（Lambert）、墨卡托投影（Mercator）、极投影。 例一： import cartopy.crs as ccrs import matplotlib.pyplot as plt #plt.axes是创建一个轴
python:Cartopy的基础使用
weixin_49284189的博客
10-30 6782
提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档 文章目录前言一、基础介绍二、使用步骤1.引入库2.读入数据总结 前言 常用地图底图的绘制一般由Basemap或者cartopy模块完成，由于Basemap库是基于python2开发的一个模块，目前已经不开发维护。故简单介绍cartopy模块的一些基础操作。 提示：以下是本篇文章正文内容，下面案例可供参考 一、基础介绍 首先导入相关模块。 import numpy as np import matplotlib.pyplot as pl
python使用cartopy画地图
weixin_41796280的博客
05-19 2073
利用cartopy进行地图绘制，利用克里金插值算法进行数据插值处理
【气象常用】wrf后处理-风场
weixin_46921265的博客
05-24 428
wrf后处理来啦
【python地理信息绘制入门】cartopy学习
qq_38734327的博客
08-18 1276
cartopy是一个用于绘制地图投影和地理数据可视化的 Python 库。它是建立在 matplotlib 的基础上，专门用于制作地图和地理数据的图表。cartopy 的目标是使地图制作变得更加简单和直观，同时提供了一些强大的功能来处理地理数据和投影。
python安装cartopy
qq_32649321的博客
07-10 895
创建新环境：https://www.machinelearningplus.com/deployment/conda-create-environment-and-everything-you-need-to-know-to-manage-conda-virtual-environment/
python cartopy绘制扇形区域图/cartopy绘制北极部分区域
weixin_43750300的博客
02-16 2107
实际上，这个问题原因还是由于投影转换的问题，在set_extend时，绘制的上下边界仍然是方形、未被正确投影的边界，与我们的set_boundary存在冲突，最根本的原因还是在于cartopy对于投影计算的一些缺陷。另外我们需要指出的是：**该方法不适用于极地投影，即NorthPolarStereo，由于NorthPolarStereo本身投影特性只需一个参数，本身并不适合。我们绘制极地投影时，同样也是使用set_boundary绘制圆形边界，那么当我们想要绘制扇形时，可以通过。请根据喜好自行选择，
cartopy：Cartopy-具有matplotlib支持的制图python库
02-25
Cartopy是一个Python软件包，旨在简化用于数据分析和可视化的绘图图。目录概述Cartopy是一个Python软件包，旨在简化用于数据分析和可视化的绘图图。 它具有以下特点： 面向对象的投影定义投影之间的点，线，多边形和...
python画图：利用matplotlib绘制函数
03-11
python画图
python基础知识笔记
03-15
python基础知识笔记
Python基础教程学习笔记 第一章 基础知识
01-05
1、python的每一个语句的后面可以添加分号也可以不添加分号；在一行有多条语句的时候，必须使用分号加以区分 2、查看Python版本号，在Dos窗口中输入“python -V”查看版本号，注意V是大写字母，这条命令是Windows...
开源的python地图
06-11
利用Python操作google map
Python基础：B站自学笔记
最新发布
07-03
跟着B站学的python基础，但是没有看完全部
python绘图 科研作图 常用函数 方法
Cherry_csc的博客
02-12 3899
python中matplotlib库中一些作图方法总结
python+matplotlib绘制带误差线的柱状图
qq_38671556的博客
01-05 2090
用python和matplotlib绘制带误差线的柱状图
小白学习cartopy气象画地图的第二天（中国区域，陆地温度分布图）
weixin_42372313的博客
02-04 9271
小白学习cartopy气象画地图的第二天（中国区域，陆地温度分布图） 首先，还是线感谢一下公众号（气象学家），因为代码和测试数据都是他的，我只是拿过来，然后慢慢读懂它，最后再以小白的角度把每一行代码解释给大家听，拾人牙慧而已。话不多说，开始！ 第一步导入库 import numpy as np import xarray as xr import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cf
python用cartopy包画地图,Python Cartopy地图，国家/地区以外的剪贴区域（多边形）...
weixin_42372837的博客
12-06 1026
I create a Stamen terrain map with the country border from NaturalEarth.Now I want to remove all the data (terrain in this case) from outside the country border.How would I do that?My example with the...
Python数据可视化——散点图、添加标签、颜色选取
小平凡的记录的博客
06-01 6247
绘制散点图、添加标签
Cartopy绘图入门指南
Jasenjane的博客
07-20 7438
嗨，你好，我是来自点点GIS的南南我与Cartopy的认识起源于“气象水文科研猫”的这个推文，那时候的我觉得，用代码画地图好酷，arcgis就感觉low爆了。但是一直没有时间学习。前段时间放暑假，磕磕绊绊装完包以后就不想动弹了，折腾环境折腾的我头皮发麻。gdal和cartopy真的是我学python以来装的最麻烦的包。很多时候我想学习，却感觉无从下手，大佬们的博文一大堆，我看了以后只有一句“卧槽，牛逼”，然后没然后了。前几天某人欺负我，把我整抑郁了，大晚上睡不着觉，就想着填个坑，虽然我技术菜，不能“炫技”，
python cartopy加图例
09-08
要在Python中使用Cartopy添加图例，你可以使用Matplotlib的colorbar()函数。下面是一个示例代码，展示了如何使用Cartopy和Matplotlib创建一个带有图例的地理热图： ``` import cartopy.crs as ccrs import matplotlib.pyplot as plt # 创建地图投影 ax = plt.axes(projection=ccrs.PlateCarree()) # 绘制地理热图 # 这里假设你已经有了一些数据，命名为data heatmap = ax.pcolormesh(lons, lats, data, transform=ccrs.PlateCarree(), cmap='hot') # 添加颜色条图例 cbar = plt.colorbar(heatmap) # 设置图例的标签 cbar.set_label('数据值') # 显示地图 plt.show() ``` 这段代码首先创建了一个地图投影对象，然后使用`pcolormesh`函数绘制地理热图。接下来，通过调用`colorbar()`函数创建一个颜色条图例对象，并使用`set_label()`方法设置图例的标签。最后，通过调用`show()`方法显示地图。 请注意，上述代码只是一个示例，你需要根据你的具体需求来适应和修改它。你可能需要调整数据的范围、颜色映射和其他绘图参数来满足你的需求。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值