通常情况下，要获取某个区域内的格点数据，如果要求不是很高，直接采取矩形框挑选方法——即锁定所需范围内的经纬度，就能挑选出需要的数据。而对于不规则的范围，数据的匹配精度有一定要求，譬如，需要严格按照某个特定区域的shapefile文件来截取数据。虽然，NCL官网提供了可行的解决方案，但是 shapefile_mask_data（包含在shapefile_utils.ncl中，官网有提供）也仅仅是较好地适用于2维的Lat-Lon数据，对于3维或者更高维度的数据，其处理效率非常低下。所以，针对于这个问题，在实际的操作中我给出了一个快速处理的方案，仅供参考：
代码语言：javascript复制 1 load "../shapefile_utils.ncl"
2
3 yearn = 14 ; 1996-2009
4 leadn = 30 ; days
5 pentadn = 6
6 itemn = 3
7 year = ispan(1996,2009,1)
8 lead = ispan(1,30,1)
9 pentad= ispan(1,6,1)
10 item = ispan(1,3,1)
11
12; Longitude co-ordinate
13; ................
14 nlon=247
15 lon=fspan(73.5,135.0,247)
16 lon!0="lon"
17 lon&lon=lon
18 lon@units="degrees_east"
19 lon@long_name="Longitude"
20 ; printVarSummary(lon)
21
22; Latitude co-ordinate
23; ................
24 nlat=169
25 lat=fspan(16.5,58.5,169)
26 lat!0="lat"
27 lat&lat=lat
28 lat@units="degrees_north"
29 lat@long_name="Latitude"
30 ; printVarSummary(lat)
31
32 filp1 = "/Users/zhpfu/Dropbox/S2S_material/code/S2S_South_China_Tibet_NEW/0.TP/"
33 f_cma = addfile(filp1 "CMA_TP_hr.nc","r")
34 f_eccc = addfile(filp1 "ECCC_TP_hr.nc","r")
35 f_ecmwf = addfile(filp1 "ECMWF_TP_hr.nc","r")
36 f_hmcr = addfile(filp1 "HMCR_TP_hr.nc","r")
37 f_ukmo = addfile(filp1 "UKMO_TP_hr.nc","r")
38
39 f_obs = addfile(filp1 "Obs_TP_hr.nc","r")
40 f_era5 = addfile(filp1 "ERA5_TP_hr.nc","r")
41 f_erai = addfile(filp1 "ERAI_TP_hr.nc","r")
42
43 tp_cma = f_cma->tp
44 tp_eccc = f_eccc->tp
45 tp_ecmwf= f_ecmwf->tp
46 tp_hmcr = f_hmcr->tp
47 tp_ukmo = f_ukmo->tp
48
49 tp_obs = f_obs->tp
50 tp_era5 = f_era5->tp
51 tp_erai = f_erai->tp
52
53 ; printVarSummary(tp_cma)
54
55;---Open shapefile and read lat/lon values.
56 dir = "/Users/zhpfu/Dropbox/S2S_material/code/South_China_Tibet/xinan/"
57 shp_filename = dir "southwest.shp"
58 ; dim_cma = dimsizes(tp_cma)
59 ; print(dim_cma)
60 ; Variable: dim_cma
61 ; Type: integer
62 ; Total Size: 20 bytes
63 ; 5 values
64 ; Number of Dimensions: 1
65 ; Dimensions and sizes: [5]
66 ; Coordinates:
67 ; (0) 27
68 ; (1) 14
69 ; (2) 30
70 ; (3) 169
71 ; (4) 247
72
73 mask_in = shapefile_mask_data(tp_erai(0,0,:,:),shp_filename,True)
74 mask_io = where(ismissing(mask_in), 0, 1)
75
76 erai_mask = tp_erai
77 erai_mask = mask (erai_mask, conform(erai_mask, mask_io, (/2,3/)), 1)
78 copy_VarCoords(tp_erai,erai_mask)
79;===================================================================
80 diro = "/Users/zhpfu/Dropbox/S2S_material/code/S2S_South_China_Tibet_NEW/0.TP/" ; Output directory
81 filo = "ERAI_TP_hr_mask.nc" ; Output file
82 system("/bin/rm -f " diro filo) ; remove if exists
83 fout = addfile (diro filo, "c") ; open output file
84 fout->tp = erai_mask其中核心代码就是：
代码语言：javascript复制1 mask_in = shapefile_mask_data(tp_erai(0,0,:,:),shp_filename,True) ; 构建母版——最基础的2D mask范围
2 mask_io = where(ismissing(mask_in), 0, 1) ;将所需范围内外数据的分离开
3
4 erai_mask = tp_erai ;简化的数组创建方案
5 erai_mask = mask (erai_mask, conform(erai_mask, mask_io, (/2,3/)), 1) ;处理高维数组进行mask
6 copy_VarCoords(tp_erai,erai_mask) ;复制坐标信息总结一下：由于使用了自带的mask、conform和where函数，相比于shapefile_mask_data基础上多层循环嵌套具有速度快、效率较高。如果你有什么更好更快的办法也欢迎留言！
—END—
