从全球数据切中国区域切割不全

星星薇

问题已经解决，使用这种内部有填充的shp文件就可以了。附件在脚本后面。
------------------------------------------------------------------------------------------------------------------------------------------
如题，使用官方的示例，从全球数据里往下切中国的nc数据，因为还要做区域平均所以这一步是必须的，脚本运行无误，但是不知道为什么相比用元数据mask的结果，切下来的数据不全，有很多地方缺失了，不知道有没有人遇到和我一样的问题。

1. ;切割nc数据切成中国区域，参考mask9.ncl,以及气象家园一些帖子。
   
2. 
   
3. load "/home/max/data/shp/china/cnmap/cnmap.ncl"
   
4. 
   
5. begin
   
6. MASK_INSIDE = False    ;  区域内还是区域外
   
7. dir = "/home/max/cmip6/fracLUT/hist/"
   
8. mpi = addfile(dir+"fracLut_Emon_MPI-ESM1-2-LR_land-hist_r1i1p1f1_gn_201001-201412.nc", "r")    ;NOTE 时间范围
   
9. mpi_fracLut = mpi ->fracLut
   
10. lat_mpi = mpi ->lat
    
11. delete(mpi)
    
12. 
    
13. ;--------- 处理missingvalue --------把nan处理掉 -------[year | 1980] x [landuse | 4] x [lat | 192] x [lon | 288]----------
    
14. if (any(isnan_ieee(mpi_fracLut))) then
    
15.       value = 1.e20
    
16.       replace_ieeenan (mpi_fracLut, value, 0)
    
17.       mpi_fracLut@_FillValue = value
    
18. end if
    
19. 
    
20. ;---读取 shp文件
    
21.   dirshp     = "/home/max/data/shp/china/cnmap/"
    
22.   f       = addfile(dirshp + "cnmap.shp", "r")
    
23.   mrb_lon = f->x
    
24.   mrb_lat = f->y
    
25.   nmrb    = dimsizes(mrb_lon)     ;计算顶点数
    
26. 
    
27.   min_mrb_lat = min(mrb_lat)    ;估计轮廓线范围
    
28.   max_mrb_lat = max(mrb_lat)
    
29.   min_mrb_lon = min(mrb_lon)
    
30.   max_mrb_lon = max(mrb_lon)
    
31. 
    
32. ;
    
33. ;实现mask效果取决于区域内外缺省值的赋值
    
34. 
    
35.   if(MASK_INSIDE) then
    
36. ;---Start with data filled in.
    
37.     mpi_fracLut_mask = mpi_fracLut   ;数据区
    
38.   else
    
39. ;---Start with data all missing
    
40.     mpi_fracLut_mask = new(dimsizes(mpi_fracLut),typeof(mpi_fracLut),mpi_fracLut@_FillValue)
    
41.     copy_VarCoords(mpi_fracLut,mpi_fracLut_mask)
    
42.   end if
    
43. 
    
44. ;让gc_inout更快的小技巧，估算轮廓线的大致范围
    
45. ;与元数据的数据范围和分辨率相同最好，速度差不了多少。  
    
46.   lat1d = fspan(-88.6, 88.6, 96)
    
47.   lon1d = fspan(0, 359, 192)
    
48.   lat1d@units = "degrees_north"
    
49.   lon1d@units = "degrees_east"
    
50. 
    
51.   ilt_mn = ind(min_mrb_lat.gt.lat1d)    ;小于轮廓线最小纬度值的索引，以下类似
    
52.   ilt_mx = ind(max_mrb_lat.lt.lat1d)  
    
53.   iln_mn = ind(min_mrb_lon.gt.lon1d)
    
54.   iln_mx = ind(max_mrb_lon.lt.lon1d)
    
55. 
    
56. ;http://bbs.06climate.com/forum.php?mod=viewthread&tid=43465&extra=&page=1
    
57. ; 已经找到错误就在下面四句，由于经纬度不同，必要时可以手动赋值。
    
58.   ilt1   = ilt_mn(dimsizes(ilt_mn)-1)    ; 纬度值起始判断点
    
59.   iln1   = iln_mn(dimsizes(iln_mn)-1)    ; 经度值起始判断点
    
60.   ilt2   = ilt_mx(0)                     ; 纬度值结束判断点
    
61.   iln2   = iln_mx(0)                     ; 经度值结束判断点
    
62. 
    
63. ;print(ilt1) ;3
    
64. ;print(ilt2) ;54
    
65. ;print(iln1) ;3
    
66. ;print(iln2) ;67
    
67. ;---蒙版区域缺省值
    
68.   if(MASK_INSIDE) then
    
69.     do ilt=ilt1,ilt2
    
70.       do iln=iln1,iln2
    
71.         if(gc_inout(lat1d(ilt),lon1d(iln),mrb_lat,mrb_lon)) then
    
72.           mpi_fracLut_mask(:,:,ilt,iln) = mpi_fracLut_mask@_FillValue
    
73.         end if
    
74.       end do
    
75.     end do
    
76.   else
    
77. ;---非蒙版区域放入数据
    
78.     do ilt=ilt1,ilt2
    
79.       do iln=iln1,iln2
    
80.         if(gc_inout(lat1d(ilt),lon1d(iln),mrb_lat,mrb_lon)) then
    
81.           mpi_fracLut_mask(:,:,ilt,iln) = mpi_fracLut(:,:,ilt,iln)
    
82.         end if
    
83.       end do
    
84.     end do
    
85.   end if
    
86. 
    
87. printVarSummary(mpi_fracLut_mask)
    
88. ; print(mpi_fracLut_mask(0,0,48,88))
    
89. ;---画图
    
90.   wks  = gsn_open_wks("eps","mask_mpi_v17")        ; send graphics to PNG file
    
91. 
    
92.   res                     = True
    
93.   res@gsnMaximize         = True           ; maximize plot in frame
    
94.   res@gsnDraw             = False          ; don't draw plot yet
    
95.   res@gsnFrame            = False          ; don't advance frame yet
    
96.   res@gsnAddCyclic        = False          ; Don't add a cyclic point.
    
97. 
    
98.   res@mpDataBaseVersion   = "MediumRes"    ; slightly better resolution
    
99. 
    
100.   res@cnFillOn=True;画填充图
     
101.   res@cnLinesOn=False;不画等值线
     
102.   res@cnLineLabelsOn=False;不要等值线上的标签
     
103. 
     
104.   res@cnFillPalette  = "cmocean_tempo"
     
105. 
     
106.   res@mpMinLatF           = 0
     
107.   res@mpMaxLatF           = 55
     
108.   res@mpMinLonF           = 70
     
109.   res@mpMaxLonF           = 140
     
110. 
     
111.   res@mpLandFillColor         = "white"
     
112. 
     
113. ;---Resources for polyline
     
114.   ; lnres                  = True
     
115.   ; lnres@gsLineColor      = "gray"
     
116.   ; lnres@gsLineThicknessF = 1.0            ; 2x thickness）
     
117. 
     
118. ; 画中国，用什么的shp画出来的就是一团乱七八糟的线
     
119.     cnres           = True
     
120.     cnres@china     = True       ;draw china map or not
     
121.     cnres@river     = False      ;draw changjiang&huanghe or not
     
122.     cnres@province  = False       ;draw province boundary or not
     
123.     cnres@nanhai    = True       ;draw nanhai or not
     
124.     cnres@diqu      = False       ; draw diqujie or not
     
125. 
     
126.    
     
127. 
     
128. ;---Draw contours of masked data
     
129.   map_mask      = gsn_csm_contour_map(wks,mpi_fracLut_mask(0,0,:,:),res)
     
130.   ;mrb_line_mask = gsn_add_polyline(wks, map_mask, mrb_lon,mrb_lat, lnres)
     
131.   chinamap = add_china_map(wks,map_mask,cnres)
     
132. 
     
133.   draw(map_mask)
     
134.   frame(wks)
     
135. 
     
136. ;----------绘制并输出第二个图
     
137.     res@tiMainString = "masked from origin data"
     
138.     res@mpOutlineOn             = True  ; Use outlines from shapefile
     
139.     res@mpDataBaseVersion       = "MediumRes"
     
140.     res@mpDataSetName           = "Earth..4"
     
141.     res@mpAreaMaskingOn         = True
     
142.     res@mpMaskAreaSpecifiers    = (/"China","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/)
     
143.     res@mpLandFillColor         = "white"
     
144.     res@mpInlandWaterFillColor  = "white"
     
145.     res@mpOceanFillColor        = "white"
     
146.     res@mpOutlineBoundarySets   = "NoBoundaries" ; 先不要画，随后再打开
     
147.     res@mpOutlineSpecifiers=(/"China","Taiwan"/)
     
148.     res@mpGeophysicalLineThicknessF=2.0 ;这两行是为了加粗边界和国界线
     
149.     res@mpNationalLineThicknessF=2.0
     
150. 
     
151. ;--------------遇到的奇怪的bug，这里的数字后面不要加.0，比如55和55.0是不同的数字，会产生错误。
     
152.     res@mpMinLatF=0
     
153.     res@mpMaxLatF=55
     
154.     res@mpMinLonF=70
     
155.     res@mpMaxLonF=140
     
156. 
     
157.     res@cnFillDrawOrder="PreDraw";先画填充
     
158.     res@cnFillOn=True;画填充图
     
159.     res@cnLinesOn=False;不画等值线
     
160.     res@cnLineLabelsOn=False;不要等值线上的标签
     
161.     res@cnFillPalette        = "cmocean_tempo"     ; set color map
     
162.   
     
163.     wks  = gsn_open_wks("eps","mask_mpi_v16")        
     
164.     plot1=gsn_csm_contour_map(wks,mpi_fracLut(0,0,:,:),res)
     
165. 
     
166.     draw(plot1)
     
167.     frame(wks)
     
168. 
     
169. 
     
170. 
     
171. 
     
172. end
     

复制代码

werewolf

为什么要切割

星星薇

werewolf 发表于 2020-11-16 17:20
为什么要切割

因为要做区域平均等处理

星星薇

俺的问题解决了，是因为shp文件的影响，换了一个shp文件就好了

某某蒙

请问您换的是哪一个shp文件呀，我也遇到了相同的问题

星星薇

换成区域内有填充的，类似这样的

星星薇

某某蒙 发表于 2020-11-21 20:03
请问您换的是哪一个shp文件呀，我也遇到了相同的问题

忘记直接回复你了QAQ，看上面

某某蒙

星星薇 发表于 2020-11-22 15:28
忘记直接回复你了QAQ，看上面

好滴好滴，我试试去，谢谢！！！

杰西卡·芬

大神你好请问也可以发我一份吗？谢谢你252277973@qq.com

星星薇

杰西卡·芬 发表于 2020-12-25 09:08
大神你好请问也可以发我一份吗？谢谢你

附件我放到后面了