shapefile,是美国环境系统研究所公司(ESRI)开发的空间数据开放格式,用于描述几何体对象:点、折线与多边形。可以存储井、河流、湖泊等空间对象的几何位置。除了几何位置,shp文件也可以存储这些空间对象的属性,例如河流的名字、城市的温度等等。
Shapefile文件指的是一种文件存储的方法,实际上该种文件格式是由多个文件组成的。其中三个文件必不可少(”.shp”, “.shx”与 “.dbf”)。表示同一数据的一组文件其文件名前缀应该相同。其中“真正”的Shapefile的后缀为shp,然而仅有这个文件数据是不完整的,必须要把其他两个附带上才能构成一组完整的地理数据。所有文件必须位于同一目录之中。
ncl可以用addfile读取shp文件,示例如下
river = True
river@gsLineThicknessF = mp_thick+1.0
river@gsLineColor = "black"
plotrv = gsn_add_shapefile_polylines(wks,plot,"/home/river.nc",river)
load "./shapefile_utils.ncl"
begin
print_shapefile_info("gadm36_CHN_2.shp") #
plot_shapefile("gadm36_CHN_2.shp") #
f = addfile("gadm36_CHN_2.shp", "r") ; Open shapefile
segments = f->segments
geometry = f->geometry
segsDims = dimsizes(segments) ;2D
geomDims = dimsizes(geometry) ;2D
geom_segIndex = f@geom_segIndex ;=0
geom_numSegs = f@geom_numSegs ;=1
segs_xyzIndex = f@segs_xyzIndex ;=0
segs_numPnts = f@segs_numPnts ;=1
lines = new(segsDims(0),graphic) ; Array to hold polygons
numFeatures = geomDims(0) ;=dimsizes(NAME_0)
lon = f->x
lat = f->y
name = f->NAME_2
segNum = 0
do i=0, numFeatures-1
startSegment = geometry(i, geom_segIndex)
numSegments = geometry(i, geom_numSegs)
print(name(i)+" startSegment:"+startSegment+" numSegments:"+numSegments)
do seg=startSegment, startSegment+numSegments-1
startPT = segments(seg, segs_xyzIndex)
endPT = startPT + segments(seg, segs_numPnts) - 1
print("startPT:"+startPT+" endPT:"+endPT)
;+" lon(startPT:endPT):"+lon(startPT:endPT)+"lat(startPT:endPT):"+lat(startPT:endPT))
segNum = segNum + 1
end do
end do
end
有一个叫regionmask的程序包,可以实现该功能。 其函数介绍页面 https://regionmask.readthedocs.io/en/stable/api.html
import cmaps, regionmask
import geopandas as gpd
ds = xr.open_dataset('~/data/t2m.nc')
var = ds['t2m'] # 三维数组
shp_file = '~/data/gadm36_CHN_0.shp'
shpf = gpd.read_file(shp_file)
mask = regionmask.mask_geopandas(shpf, ds)
var = var.where(mask==0.0).mean(('lat','lon')) # 虽然var和mask的维数不同,但不影响
# ds也可以换成var,只要是包含经纬度信息的xarray数据就可以
# 且默认ds中经纬度名字为lat,lon,若不是,需要指定,此时语法是
mask = regionmask.mask_geopandas(shpf, ds, 'longitude', 'latitude')
# 也可以在shpf后直接跟lon和lat数组,且该数组可以是一维或二维,可以是xarray或array
mask = regionmask.mask_geopandas(shpf, ds['longitude'], ds['latitude'])
# 无论哪种方法,都只返回一个二维数组xarray.DataArray mask,其中不位于shp_file定义的区域内的格点为nan
# 由于此处shp_file只定义了一个区域,所以在区域中的格点为0
# 若shp_file有多个区域,则mask中的值从0-len(区域)
# regionmask中已包含了许多常见的区域划分,可直接调用,不需要使用shapefile
# natural_earth_v5_0_0.countries_110中包含177个区域,其中China的编号是139
region = regionmask.defined_regions.natural_earth_v5_0_0.countries_110
mask = region.mask(var) #输入经纬度信息的方法和返回数组信息同regionmask.mask_geopandas
var = var.where(mask==139)
nan: not a number
inf: 无穷大
xarray的DataArray的求和求平均函数会自动忽略 nan。但若有 inf,则求和与求平均的结果都是 inf。
numpy中大多数函数都不会忽略 nan 和 inf。若数组中有 nan 和 inf, 求和求平均的结果都是 nan or inf。举例如下
a = np.array([1,2,3,np.nan,4,5])
print('%f, %f'%(a.min(),a.max())) # 结果为 nan,nan
print('%f, %f'%(np.nanmin(a),np.nanmax(a))) # 结果为1,5
print('%f, %f'%(a.sum(),a.mean())) # 结果也是 nan,nan
print('%f, %f'%(np.nansum(a),np.nanmean(a))) # 结果为 15, 3
a = np.array([1,2,np.nan,np.inf,4,5])
print('%f, %f'%(a.min(),a.max())) # 结果为 nan,nan
print('%f, %f'%(np.nanmin(a),np.nanmax(a))) # 结果 1, inf
print('%f, %f'%(a.sum(),a.mean())) # 结果 nan,nan
print('%f, %f'%(np.nansum(a),np.nanmean(a))) # 结果 inf, inf
a = np.array([1,2,3,6,4,5])
a = np.ma.array(a,mask=(a==5)) # 结果 [1.0 2.0 3.0 6.0 4.0 --]
print('%f, %f'%(a.sum(),a.mean())) # 结果 16.000000, 3.200000
print('%f, %f'%(np.nansum(a),np.nanmean(a))) # 结果同上
# 使用scipy.stats.计算相关系数时,不能有缺测存在,因此需要替换缺测
# 可以自己指定将nan,np.inf替换为某个具体数值,或者使用默认的数值
new = np.nan_to_num(a, nan=0, posinf=999, neginf=-999)
np.sqrt(a), 若a是实数且a中有负数,则负数的计算值为nan,若a是复数且有负数,则负数的计算结果为复数
a = np.array([1,2,np.nan,np.inf,4,5])
b = np.ma.array(a,mask=(a==4))
# 得到的b是
# masked_array(data=[1.0, 2.0, nan, inf, --, 5.0],
# mask=[False, False, False, False, True, False],
# fill_value=1e+20)
da = xr.DataArray(np.arange(16).reshape(4, 4), dims=["x", "y"])
da.where(da.x + da.y < 4) # 保留满足条件的变量,未满足条件的变量设为nan
# da.where(da.x + da.y < 4,0),会将未满足条件的变量设为0
# output is
# <xarray.DataArray (x: 4, y: 4)>
# array([[ 0., 1., 2., 3.],
# [ 4., 5., 6., nan],
# [ 8., 9., nan, nan],
# [12., nan, nan, nan]])
# Dimensions without coordinates: x, y
最近发现 python 一个非常操蛋的内存存储问题,好像叫 reference count,一不留心就容易计算失误。
import numpy as np
def diff(a):
b=a
b[1:4]=a[1:4]-a[2:]
return b
a=np.array([0,2,1,3,4])
c=diff(a)
print(a) # array([ 0, 1, -2, -1, 4])
print(c) # array([ 0, 1, -2, -1, 4])
可以看到变量a也发生了变化。问题在于b=a
这个命令在python中不属于赋值,只是为a多增加了一个索引,a和b都指向同一个内存地址。因此改变b的同时,a也会发生变化。
b=a[0:2]
numpy数组的切片赋值也存在类似的问题。但列表的切片赋值不存在这个问题。
避免这个问题的解决方法是b=a.copy()
。或者a进行简单运算后再赋值给b,此时a与b指向不同的内存地址,如下例所示。
b=a/1.0 #
b[0]=2
print(b) # array([ 2., 1., -2., -1., 4.])
print(a) # array([ 0, 1, -2, -1, 4])
xarray中也存在类似的问题,解决方法b=foo[1:4].data.copy()
。
import xarray as xr
a=np.array([0,2,1,3,4])
foo = xr.DataArray(a, dims=["time"])
b=foo[1:4].data
b[0]=10
print(foo) # foo也发生了变化
# <xarray.DataArray (time: 5)>
# array([ 0, 10, 1, 3, 4])
# Dimensions without coordinates: time
print(b)
# [10 1 3]
numpy主要用于数组与矩阵运算,主要数据结构就是N维同类型数组array,可以用array.shape
, array.dtype
查看数组结构和类型
a = np.array([1,2,3],dtype=float)
b = np.zeros(a.shape)
c = np.zeros( [len(lev),len(behv),len(months),len(year)],dtype=float )
x = np.arange(6).reshape(2,3) + 10
# 得到数组 array([[10, 11, 12],[13, 14, 15]])
ind = np.argwhere(x>11)
# 得到一个二维数组,表示满足条件的位置坐标
# array([[0, 2],[1, 0],[1, 1],[1, 2]])
# 此时 x[ind[:,0],ind[:,1]] 为 array([12, 13, 14, 15])
ind = np.unravel_index(np.argmax(x, axis=None), x.shape)
# 得到的ind 为(1,2), x[ind] = 15
# 若直接 np.argmax(x, axis=None),得到的结果是5
# 若 np.argmax(x, axis=0),得到的结果 array([1, 1, 1])
两个不同维度的数组(a/b)相除,只有当b的维数与a的最后几个维数相等时才可以相除
a=np.arange(24).reshape((4, 3, 2))
b=np.arange(6).reshape((3, 2))
print(a/b)
#array([[[ nan, 1. ],[ 1. , 1. ],[ 1. , 1. ]],
# [[ inf, 7. ],[ 4. , 3. ],[ 2.5, 2.2]],
# [[ inf, 13. ],[ 7. , 5. ],[ 4. , 3.4]],
# [[ inf, 19. ],[10. , 7. ],[ 5.5, 4.6]]])
#
b=np.array([1,2,3,4])
c=a/b # 报错
# ValueError: operands could not be broadcast together with shapes (4,3,2) (4,)
# 此时可以将a中的第一维度移动到最后,这样就可以和b相除
c = np.moveaxis(a,0,-1)/b # c.shape 等于 (3, 2, 4)
d = np.moveaxis(c,-1,0) # d.shape 等于 (4, 3, 2)
# 也可以将b变成和a一样维数的数组,然后再相除,得到的结果一样
d = a/np.broadcast_to(b.reshape(4, 1, 1), a.shape)
a=np.array([tim3.time.dt.year.isin(1990),tim3.time.dt.month.isin(2)])
print(a)
print(a.all(axis=0))
# output:
#[[ True True False]
# [False True False]]
#
#[False True False]
new = np.where(a<5, a, 10*a)
# 原数组为 array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
# 新数组为 array([0, 1, 2, 3, 4, 50, 60, 70, 80, 90])])
a = xr.DataArray(np.arange(25).reshape(5, 5), dims=("x", "y"))
new = xr.where(a<5, a, 10*a) # 用法通np.where
# 但 xarray.DataArray.where 的用法和 xr.where 不同
a.where(a.x + a.y < 5, -1) # 将不满足条件的值设为-1,若不给定,则会设为缺测
a.where(a.x + a.y < 5, drop=True) # 将不满足条件的值设为缺测,并删去都是缺测的列或行
# a.x + a.y < 5返回逻辑数组,因此也可以用 a.x.isin([2,3]) 代替
# numpy 也有 isin的用法,但格式和 xarray 不一样
b = np.array([range(5,15,1),range(0,10,1)])
mask = np.isin(b, range(5,10))
# mask是一个shape同b的逻辑数组,
# b[mask]可以得到 array([5, 6, 7, 8, 9, 5, 6, 7, 8, 9])
mask_invert = np.isin(b, range(5,10), invert=True)
# mask_invert是mask的逻辑反面,即mask为True的地方,mask_invert为False
numpy.savetxt(fname,X,fmt='%.18e')
# 第一个参数为文件名,第二个参数为需要存的数组(一维或者二维)
# 数组的第一个维度变为行,第二个维度变为列
# 这个写入功能用得比较少,感觉没有python自带的写入功能好用
# 详情可参考 https://renqlsysu.github.io/2021/10/04/python_file/
data = numpy.loadtxt(fname,dtype=<class 'float'>,
comments='#', delimiter=None, usecols = (0,1))
# 将数据读出为array类型,若文件有ny行,nx列,得到的数组维度是(ny,nx)
# 可以用skiprows跳过前几行,用usecols指定读取的列(从0列开始)
Matplotlib有很多自带的色标,但感觉都比较丑, 作为用惯ncl的人,还是觉得ncl的色标比较舒适,因此网上有大神写了一个ncl色标的python库cmaps。
from matplotlib import cm
import cmaps
ncmap = cm.viridis # <matplotlib.colors.ListedColormap object at 0x2b6b1578b550>
ncmap = cm.get_cmap('viridis',17) # 从viridis中取出17个颜色,此时ncmap.N是17
ncolor = cm.viridis.N #会现实该色标的颜色数量,256
ncmap = cmaps.precip3_16lev # <cmaps.colormap.Colormap object at 0x2b6b0abdceb0>,
# 也可以用ncmap.N统计其颜色数量,该色标有17个颜色
# cmaps.BlueDarkRed18.N的颜色有18个颜色
#可以通过以下方法获得colormap的RGBA值,一个Nx4的array数组,N是该colormap的颜色数量,4代表各个颜色的RGBA值
data = ncmap.colors
data = ncmap(range(8)) #取前8个颜色
data = ncmap(np.linspace(0, 1, 8)) #显示线性插值的8个颜色
from matplotlib.colors import ListedColormap
ncmap = ListedColormap(["blue","green","red"])
# 合并两个色标
newcolors = np.vstack((cm.Oranges_r(np.linspace(0, 1, 128)),
cm.Blues(np.linspace(0, 1, 128))))
newcmp = ListedColormap(newcolors, name='OrangeBlue')
# 将已定义的色标转置
ncmap = ListedColormap(cmaps.precip3_16lev(range(0,17,1))[::-1])
# 获得较为稀疏的已定义坐标
ncmap = cm.get_cmap('viridis',17) # matplotlib自带的colormap可以这样操作
ncmap = ListedColormap(cmaps.GMT_globe(range(0,209,3)))
# 替换已定义色标中的某一颜色
newcolors = cm.viridis(np.linspace(0, 1, 256))
pink = np.array([248/256, 24/256, 148/256, 1])
newcolors[:25, :] = pink
newcmp = ListedColormap(newcolors)
# 将值线性分布到colormap上
pcm = ax.pcolormesh(x, y, Z, vmin=-1., vmax=1., cmap='RdBu_r')
# 将值对数分布到colormap上
pcm = ax[0].pcolor(X, Y, Z,
norm=colors.LogNorm(vmin=Z.min(), vmax=Z.max()),
cmap=cmaps.BlueDarkRed18, shading='auto')
# 将值线性分布到colormap上,同时vcenter的值将位于colormap中间
pc = ax2.pcolormesh(Z, norm=colors.CenteredNorm(vcenter=0), cmap=cm.coolwarm)
# 定义离散的colormap值
# 一般我会将ncolors设为len(bounds)+1,但看这个例子好像ncolors只要大于len(bounds)就可以
import matplotlib.colors as colors
bounds = np.array([-0.2, -0.1, 0, 0.5, 1])
norm = colors.BoundaryNorm(boundaries=bounds, ncolors=cm.viridis.N, extend='both')
pcm = ax[2].pcolormesh(X, Y, Z, norm=norm, cmap=cm.viridis)
# 使值在colormap中间两侧的值的变化率不同
# make a colormap that has land and ocean clearly delineated and of the
# same length (256 + 256)
colors_undersea = plt.cm.terrain(np.linspace(0, 0.17, 256))
colors_land = plt.cm.terrain(np.linspace(0.25, 1, 256))
all_colors = np.vstack((colors_undersea, colors_land))
terrain_map = colors.LinearSegmentedColormap.from_list(
'terrain_map', all_colors)
# make the norm: Note the center is offset so that the land has more
# dynamic range:
divnorm = colors.TwoSlopeNorm(vmin=-500., vcenter=0, vmax=4000)
pcm = ax.pcolormesh(longitude, latitude, topo, rasterized=True, norm=divnorm,
cmap=terrain_map, shading='auto')
cb = fig.colorbar(pcm, shrink=0.6)
cb.set_ticks([-500, 0, 1000, 2000, 3000, 4000])
参考:
目前国际上使用的大气再分析资料已经历了多次迭代,除了以往被人们所熟知的NCEP1、NCEP2、ERA-interim,现在又出现了以下几个高分辨率的新品种(按发布顺序)。
NCEP(National Center for Environmental Prediction)于2010年发布的第3代全球再分析资料,全称Climate Forecast System Reanalysis,使用Coupled Forecast System (CFS) model,a coupled atmosphere–ocean–sea ice–land model, assimilates conventional observations and satellite radiance, and includes the time evolution of CO2 concentrations. 空间分辨率T382,约38km,垂直64层,使用GSI 3DVAR同化方法。
逐6h数据下载:https://rda.ucar.edu/datasets/ds093.0/#!description
日本气象厅(Japanese Meteorological Agency)2013年发布的第二代再分析数据,基于全球谱模式JMA2002,水平分辨率T319(差不多0.5°分辨率),垂直60层,4D-Var semi-Lagrangian assimilation schemes with Variational Bias Correction (VarBC) for satellite radiances,采用新的辐射方案,引入随时间变化的温室气体浓度。有逐3h,逐6h。
数据下载地址:https://rda.ucar.edu/datasets/ds628.0/
日本气象厅的官网介绍:https://jra.kishou.go.jp/JRA-55/index_en.html
Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA-2)由 National Aeronautics and Space Administration (NASA) Global Modeling and Assimilation Office (GMAO)于2015年发布,基于Goddard Earth Observing System Model, Version 5 (GEOS-5) data assimilation system,时间跨度为1980至现在,分辨率 0.5lat x 0.625lon,垂直72层,有逐小时,逐3h,逐6h。
存在的缺陷:极低海域降水太大,热带区域山地降水偏多
数据官网介绍: https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/data_access/
欧洲中期天气预报中心(European Centre for Medium-range Weather Forecasts,ECMWF)于2020年发布的、逐小时、31km水平分辨率的第五代再分析资料
我国国家气象信息中心于2021年5月发布的、逐6小时、 34km水平分辨率的我国第一代全球大气和陆面再分析产品(CRA),下载网址 https://data.cma.cn/data/cdcdetail/dataCode/NAFP_CRA40_FTM_6HOR.html,但感觉下载不是很方便,需要注册一个国家气象科学数据中心的教育科研实名注册用户(需要上传省部级以上的科研项目书)。但气象局内部的人下载可能会方便一些。
我国的该项目启动于2014年,历时6年研制。基于Global System Model (GSM) of the Global Forecast System (GFS) of NCEP and its Gridpoint Statistical Interpolation (GSI) 3DVAR data assimilation system,同化了多种传统观测和卫星观测数据,尤其是东亚。
根据使用过的人反馈,质量比NCEP1好。但和ERA5及JRA55没有做过详细对比。
Li, C., T. Zhao, C. Shi, and Z. Liu, 2021: Assessment of precipitation from the CRA40 dataset and new generation reanalysis datasets in the global domain. Int. J. Climatol., 41, 5243–5263.
Yu, X., L. Zhang, T. Zhou, and J. Liu, 2021: The Asian subtropical westerly jet stream in CRA-40, ERA5, and CFSR reanalysis data: Comparative assessment. J. Meteorol. Res., 35, 46–63.