ncl 添加点shp文件_气象编程 | NCL高效快速精准提取不规则区域内的格点数据

添加新云天气象主编微信或QQ：130188121，及时获取或发布气象升学、就业、会议、征稿及学术动态等信息！

通常情况下，要获取某个区域内的格点数据，如果要求不是很高，直接采取矩形框挑选方法——即锁定所需范围内的经纬度，就能挑选出需要的数据。而对于不规则的范围，数据的匹配精度有一定要求，譬如，需要严格按照某个特定区域的shapefile文件来截取数据。虽然，NCL官网提供了可行的解决方案，但是 shapefile_mask_data(包含在shapefile_utils.ncl中，官网有提供)也仅仅是较好地适用于2维的Lat-Lon数据，对于3维或者更高维度的数据，其处理效率非常低下。所以，针对于这个问题，在实际的操作中我给出了一个快速处理的方案，仅供参考：

 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

其中核心代码就是：

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)  ;处理高维数组进行mask6   copy_VarCoords(tp_erai,erai_mask)                 ;复制坐标信息

总结一下：由于使用了自带的mask、conform和where函数，相比于shapefile_mask_data基础上多层循环嵌套具有速度快、效率较高。如果你有什么更好更快的办法也欢迎留言！

有任何问题都欢迎交流探讨，共同学习进步！