目录
【总代码】(简易好理解版本,最后面还有个小升级版)
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import netCDF4 as nc
from pyproj import Proj
import datetime as dt
import xarray as xr
import datetime as dt
class Sim1(object):
def __init__(self,file_nc):
self.file_obj = xr.open_dataset(file_nc)
pass
def Draw_vector(self,mon,n):
u_mean,v_mean,speed = self.cal_mean(mon,n)
m,xi,yi = self.bmap()
m.quiver(xi[::4,::4],yi[::4,::4],u_mean[::4,::4],v_mean[::4,::4],speed[::4,::4],scale=180,color='r',width=0.002)#scale=80--80m/s设置为箭头单位强度
plt.savefig(f'sim_{n}.png')
pass
def cal_mean(self,mon,n):
u = self.file_obj.u
v = self.file_obj.v
u = u.loc[u.time.dt.month.isin(mon)]
v = v.loc[v.time.dt.month.isin(mon)]
u_mean = np.mean(u,axis=0)
v_mean = np.mean(v,axis=0)
speed = np.sqrt(u_mean**2 + v_mean**2)
return u_mean,v_mean,speed
def bmap(self):
m = Basemap(boundinglat=75, lon_0 = 0,projection = 'npstere', round=True)
m.drawcoastlines(linewidth=0.5)
m.drawstates()
m.drawcountries()
lons = self.file_obj.longitude
lats = self.file_obj.latitude
xi,yi = m(lons,lats)
return m,xi,yi
a = Sim1('/mnt/g/st-nsidc/icemotion_daily_nh_25km_20210101_20211231_v4.1.nc')
a.Draw_vector(mon=range(1,13),n='2021')
【解析】
一、相关库的安装
首先大家的环境要有这些库
可以用conda list查询:
如果没有(例如没有netCDF4)则用pip install netCDF4或者conda install netCDF4安装对应的库。
二、数据下载
进入NSIDC官网:Homepage | National Snow and Ice Data Center (nsidc.org)
点击Explore Data
输入要下载的产品的名字:如sea ice motion
点击选择NetCDF格式,选择一个Parameter是SEA ICE MOTION的产品即可
(我选择了第一个Polar Pathfinder)
点进去往下滑,点击HTTPS
有南极北极两个文件夹,这里选择北极north
里面有四个文件夹,主要用的是daily(每日数据)和weekly(每周数据)
点击要下载的文件就好啦
三、代码解析
#定义一个类,这个类的功能就是完成海冰运动矢量制图
class Sim1(object):
def __init__(self,file_nc):#初始化类,传入文件路径
self.file_obj = xr.open_dataset(file_nc)#将文件传给对象的file_obj这个属性
pass
def Draw_vector(self,mon,n):#1️⃣画矢量图函数
u_mean,v_mean,speed = self.cal_mean(mon,n) #算出U,V的平均值和合成速度(具体函数看下面)
m,xi,yi = self.bmap() #得到地理地图m,和每个u,v经过投影后的xi,yi(具体函数看下面)
m.quiver(xi[::4,::4],yi[::4,::4],u_mean[::4,::4],v_mean[::4,::4],speed[::4,::4],scale=180,color='r',width=0.002)#quiver函数画图,xi,yi为点的坐标,u_mean,v_mean为平均值,speed是合成速度,scale是设置箭头单位强度,color是设置颜色,width时设置箭头宽度,更多参数......
plt.savefig(f'sim_{n}.png') #Linux系统用这个,Win用plt.show()
pass
def cal_mean(self,mon,n):2️⃣计算平均值函数
u = self.file_obj.u
v = self.file_obj.v
u = u.loc[u.time.dt.month.isin(mon)] #将月份为mon的U选出来
v = v.loc[v.time.dt.month.isin(mon)]
u_mean = np.mean(u,axis=0) #计算选出来的u的平均值,axis=0是对第一维度求平均,axis=1对第二维度求平均,这里u的时间对应是第一维,所以是对第一维(时间)求平均
v_mean = np.mean(v,axis=0)
speed = np.sqrt(u_mean**2 + v_mean**2) #勾股定理求速度值的大小
return u_mean,v_mean,speed
def bmap(self):3️⃣创建地图
m = Basemap(boundinglat=75, lon_0 = 0,projection = 'npstere', round=True) #创建地理底图
m.drawcoastlines(linewidth=0.5)#这三行是添加相应的地理要素 海洋线、国界等
m.drawstates()
m.drawcountries()
lons = self.file_obj.longitude #将每个点的经纬度取出
lats = self.file_obj.latitude
xi,yi = m(lons,lats) #这样写就能转成m投影下对应的坐标了
return m,xi,yi
a = Sim1('/mnt/g/st-nsidc/icemotion_daily_nh_25km_20210101_20211231_v4.1.nc') #文件路径
a.Draw_vector(mon=range(1,13),n='2021') #mon输入要求平均的月份,n是图像名字的后缀(如mon=345,n='spring')
运行结果如下:
年平均(mon=rang(1,13),n='2021')
春季平均(mon=(3,4,5),n='spring')
四、代码改进(小升级版本)
功能实现:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import netCDF4 as nc
from pyproj import Proj
import datetime as dt
import xarray as xr
import datetime as dt
import cartopy.crs as ccrs
class Sim1(object):
def __init__(self,file_nc):
self.file_obj = xr.open_dataset(file_nc)
pass
def Draw_vector(self,mon,n):
u_mean,v_mean,speed = self.cal_mean(mon,n)
m,xi,yi = self.bmap()
m.quiver(xi[::4,::4],yi[::4,::4],u_mean[::4,::4],v_mean[::4,::4],speed[::4,::4],scale=180,color='r',width=0.002)#scale=80--80m/s设置为箭头单位强度
plt.title(f'{n}',x=0.5,y=-0.15 ,color='b',size=7.5)
#plt.savefig(f'sim_{n}.png')
pass
def cal_mean(self,mon,n):
u = self.file_obj.u
v = self.file_obj.v
u = u.loc[u.time.dt.month.isin(mon)]
v = v.loc[v.time.dt.month.isin(mon)]
u_mean = np.mean(u,axis=0)
v_mean = np.mean(v,axis=0)
speed = np.sqrt(u_mean**2 + v_mean**2)
return u_mean,v_mean,speed
def bmap(self):
m = Basemap(boundinglat=75, lon_0 = 0,projection = 'npstere', round=True)
m.drawcoastlines(linewidth=0.5)
m.drawstates()
m.drawcountries()
m.drawmeridians(np.arange(0, 360, 30),labels=[True,True,False,False],size=5, latmax=90)
m.drawparallels(np.arange(75, 90, 5),labels=[False,False,True,False],size=5)
lons = self.file_obj.longitude
lats = self.file_obj.latitude
xi,yi = m(lons,lats)
return m,xi,yi
调用函数:
a = Sim1('/mnt/g/st_NSIDC/icemotion_daily_nh_25km_20210101_20211231_v4.1.nc')
fig = plt.figure(figsize=[6, 6])
ax1 = fig.add_subplot(2, 2, 1)
ax1 = a.Draw_vector(mon=(3,4,5),n='spring')
ax2 = fig.add_subplot(2, 2, 2)
ax2 = a.Draw_vector(mon=(6,7,8),n='summer')
ax3 = fig.add_subplot(2, 2, 3)
ax3 = a.Draw_vector(mon=(9,10,11),n='autumn')
ax4 = fig.add_subplot(2, 2, 4)
ax4 = a.Draw_vector(mon=(1,2,12),n='winter')
fig = fig.tight_layout(pad=2, w_pad=0.8, h_pad=0.8)
plt.suptitle('SIM_2021')
##主图titlell???
plt.savefig('SIM_2021.png',dpi=400)
结果图: