Python使用cartopy库画南极站点的等高线热力图

       最近课题需要画南极科考站的热力图，因为之前都用MATLAB画图，没怎么用过Python，从安装cartopy、到画地图、画站点、画热力图走了很多的坑，所以整理了一下从安装库到画热力图的过程。

1.安装cartopy库到虚拟环境

       安装cartopy库这里卡了好久，最开始是用pip安装，发现pip安装以后还自己需要安装cartopy的依赖库，而且安装的依赖库版本也有限制，最麻烦的是安装完以后还与自己之前安装的库有冲突，后来看官方文档，才知道用conda可以直接安装。为了避免与之前安装的库有冲突，所以最好新建一个虚拟环境，再使用conda命令安装。本人使用的Jupyter笔记本，使用其他IDE原理基本差不多，但使用Jupyter会方便很多。

1.打开Anaconda Prompt，用conda创建虚拟环境，可指定Python版本：

conda create -n myenv python=3.9

2.进入创建的虚拟环境：

activate myenv

3.安装ipykernel包：

pip install --user ipykernel

4.将虚拟环境加入Jupyter：

python -m ipykernel install --user --name=myenv

5.运行jupyter notebook或jupyter lab，即可使用刚刚创建的虚拟环境。

2.在虚拟环境中安装cartopy库

     建议使用官网推荐的conda命令会简单很多：

conda install -c conda-forge cartopy

     输入命令以后会自动安装cartopy库。

     安装完成以后，输入如下代码，如果安装成功将会出现matplotlib画的世界地图

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

ax = plt.axes(projection=ccrs.Mollweide())
ax.stock_img()
plt.show()

3.画南极地图

fig = plt.figure(figsize=[13, 7])#确定图片的大小
ax = fig.add_subplot(1, 1, 1, projection=ccrs.SouthPolarStereo())
ax= plt.axes(projection=ccrs.SouthPolarStereo())#选择投影方式
ax.set_extent([-180, 180, -90, -60], ccrs.PlateCarree()) # 设置显示经度、纬度的范围
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.gridlines()

plt.show()

       得到SouthPolarStereo()投影方式的南极地图。

 4.根据科考站的经纬度画等高线热力图

        科考站的位置信息可以在百度百科上查到，画热力图的方法：首先在画布上生成均匀的网格，再根据每一个具体站点的位置坐标使用二维概率密度函数生成以该站点位置坐标为均值的高斯概率密度函数，最后在将所有的概率密度函数的值相加，即可画出所有站点的等高线热力图。这里要注意一点，站点的坐标最好不要距离太近，太近的话会导致该点的概率密度值特别大，画热力图的时候导致其他的点会变得很淡。

在网格上生成以坐标点为均值的二维高斯下降的概率密度函数。

from scipy.stats import multivariate_normal

M = 3600  #3600x3600的网格
x = np.linspace(-180, 180, 3600) #经度
y = np.linspace(-90,-30,3600) #维度
X,Y = np.meshgrid(x,y) #将x和y网格化
d = np.dstack([X,Y])

z_1=multivariate_normal(mean=(77.116449,-80.41734), cov=[[10,0],[0,0.5]])#Kunlun Station
z_2=multivariate_normal(mean=(76.966667,-73.85), cov=[[10,0],[0,0.5]]) #Taishan Station
z_3=multivariate_normal(mean=(76.371652,-69.373587), cov=[[10,0],[0,0.5]])#Zhongshan

Z1 = z_1.pdf(d).reshape(M,M)
Z2 = z_2.pdf(d).reshape(M,M)
Z3 = z_3.pdf(d).reshape(M,M)

Z =  Z1 + Z2 + Z3

    画等高线热力图

fig = plt.figure(figsize=[13, 7])
ax = fig.add_subplot(1, 1, 1, projection=ccrs.SouthPolarStereo())
ax= plt.axes(projection=ccrs.SouthPolarStereo())
ax.set_extent([-180, 180, -90, -65], ccrs.PlateCarree()) # 经度 纬度 范围


filled_c = ax.contourf(X, Y, Z, cmap='Spectral_r',transform=ccrs.PlateCarree(), alpha = 0.7)
#line_c = ax.contour(X, Y, Z, levels=filled_c.levels,colors=['blue'],transform=ccrs.PlateCarree())

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE,lw=0.25)
ax.gridlines()

画出的南极站点热力图如下。

完整代码

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from scipy.stats import multivariate_normal

M = 3600
x = np.linspace(-180, 180, 3600)
y = np.linspace(-90,-30,3600)
X,Y = np.meshgrid(x,y)
d = np.dstack([X,Y])

z_1=multivariate_normal(mean=(77.116449,-80.41734), cov=[[10,0],[0,0.5]])#Kunlun
z_2=multivariate_normal(mean=(-34.627588,-77.873696), cov=[[10,0],[0,0.5]])#Belgrano II
z_3=multivariate_normal(mean=(-26.197713,-75.612501), cov=[[10,0],[0,0.5]])#Halley
z_4=multivariate_normal(mean=(39.581836,-69.004122), cov=[[10,0],[0,0.5]])#Showa
z_5=multivariate_normal(mean=(62.873726,-67.602746), cov=[[10,0],[0,0.5]])#Mawson
z_6=multivariate_normal(mean=(76.966667,-73.85), cov=[[10,0],[0,0.5]]) #Taishan Station
z_7=multivariate_normal(mean=(2.535138,-72.011662), cov=[[10,0],[0,0.5]])#Troll
z_8=multivariate_normal(mean=(11.731944,-70.766667), cov=[[10,0],[0,0.5]])#Maitri
z_9=multivariate_normal(mean=(93.009724,-66.553122), cov=[[10,0],[0,0.5]])#Mirny
z_10=multivariate_normal(mean=(77.9675,-68.576667), cov=[[10,0],[0,0.5]])#Princess Elizabeth Land
z_11=multivariate_normal(mean=(0.066792,-75.001882), cov=[[10,0],[0,0.5]])#Kohnen
z_12=multivariate_normal(mean=(164.228815,-74.624015), cov=[[10,0],[0,0.5]])#Jang Bogo
z_13=multivariate_normal(mean=(110.526613,-66.282514), cov=[[10,0],[0,0.5]])#Casey
z_14=multivariate_normal(mean=(140.001111,-66.662778), cov=[[10,0],[0,0.5]])#Dumont d'Urville Station
z_15=multivariate_normal(mean=(123.332196,-75.09978), cov=[[10,0],[0,0.5]])#Concordia
z_16=multivariate_normal(mean=(166.668235,-77.846323), cov=[[10,0],[0,0.5]])#McMurdo
z_17=multivariate_normal(mean=(-44.737891,-60.737963), cov=[[10,0],[0,0.5]])#Orcadas
z_18=multivariate_normal(mean=(-83.261666,-79.768036), cov=[[10,0],[0,0.5]])#Union Glacier
z_19=multivariate_normal(mean=(106.837328,-78.464422), cov=[[10,0],[0,0.5]])#Vostok

Z1 = z_1.pdf(d).reshape(M,M)
Z2 = z_2.pdf(d).reshape(M,M)
Z3 = z_3.pdf(d).reshape(M,M)
Z4 = z_4.pdf(d).reshape(M,M)
Z5 = z_5.pdf(d).reshape(M,M)
Z6 = z_6.pdf(d).reshape(M,M)
Z7 = z_7.pdf(d).reshape(M,M)
Z8 = z_8.pdf(d).reshape(M,M)
Z9 = z_9.pdf(d).reshape(M,M)
Z10 = z_10.pdf(d).reshape(M,M)
Z11 = z_11.pdf(d).reshape(M,M)
Z12 = z_12.pdf(d).reshape(M,M)
Z13 = z_13.pdf(d).reshape(M,M)
Z14 = z_14.pdf(d).reshape(M,M)
Z15 = z_15.pdf(d).reshape(M,M)
Z16 = z_16.pdf(d).reshape(M,M)
Z17 = z_17.pdf(d).reshape(M,M)
Z18 = z_18.pdf(d).reshape(M,M)
Z19 = z_19.pdf(d).reshape(M,M)


Z =  Z1 +Z2 + Z3 + Z4+ Z5  +Z6 +Z7 +Z8 + Z9+ Z10 +Z11 + Z12+ Z13 + Z14 + Z15 + Z16 + Z17 + Z18 +Z19

fig = plt.figure(figsize=[13, 7])
ax = fig.add_subplot(1, 1, 1, projection=ccrs.SouthPolarStereo())
ax= plt.axes(projection=ccrs.SouthPolarStereo())
ax.set_extent([-180, 180, -90, -65], ccrs.PlateCarree()) # 经度 纬度 范围


filled_c = ax.contourf(X, Y, Z, cmap='Spectral_r',transform=ccrs.PlateCarree(), alpha = 0.5)
#line_c = ax.contour(X, Y, Z, levels=filled_c.levels,colors=['blue'],transform=ccrs.PlateCarree())

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE,lw=0.25)
ax.gridlines()