NCL 站点作图插值

546243727

四川省156个站的暴雨日数空间分布图、脚本、数据如是
请教各位大神：
1.  olon = fspan(97,109,130） ，olat = fspan(26,35,100)  控制了分辨率吗 ， 是把分辨率设置成0.1x0.1了吗？

2.  数据是四川156个站的站点数据，为什么插出省外如此多？ 是 rscan = (/10,7,4,1/)  和
final = obj_anal_ic_deprecated(lon,lat,baoyu,olon,olat,rscan,False)  用得不对吗？

3.  数据没有缺测是不是不用写baoyu@_FillValue = -999.000000？

4. 怎么MARK掉四川以外的地区呢？  

谢谢各位大神！


load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"
begin
  data  = asciiread("/home/lhz/sta_pre_days/baoyu.txt",156,"float")
  sta   = asciiread("/home/lhz/sta_pre_days/station156.txt",(/156,3/),"float")
  baoyu = data(:)
  lon   = sta(:,1)
  lat   = sta(:,2)
  olon = fspan(97,109,130)
  olat = fspan(26,35,100)

  olon!0          = "lon"
  olon@long_name  = "lon"
  olon@units      = "degrees-east"
  olon&lon        = olon
  olat!0          = "lat"
  olat@long_name  = "lat"
  olat@units      = "degrees_north"
  olat&lat        = olat

  baoyu@_FillValue = -999.000000
  rscan = (/10,7,4,1/)   
  final = obj_anal_ic_deprecated(lon,lat,baoyu,olon,olat,rscan,False)

  wks = gsn_open_wks ("png","baoyu")  
   
  res                         = True            
  res@gsnMaximize             = True
  res@gsnDraw                 = False
  res@gsnFrame                = False

;>--------------------------------------------<
;            set for the map
;>--------------------------------------------<
  res@mpMinLatF               = 26.                        
  res@mpMaxLatF               = 35.
  res@mpMinLonF               = 97.
  res@mpMaxLonF               = 109.
  res@tmXBMode                = "Explicit"   
  res@tmXBValues              = (/97,99,101,103,105,107,109/)
  res@tmXBLabels              = (/"97~S~o~N~E","99~S~o~N~E","101~S~o~N~E","103~S~o~N~E","105~S~o~N~E","107~S~o~N~E","109~S~o~N~E"/)   
  res@tmXBMinorValues         = fspan(97,109,31)
  res@tmXBMinorOn             = True  
  
  res@tmYLMode                = "Explicit"   
  res@tmYLValues              = (/26,28,30,32,34,36/)
  res@tmYLLabels              = (/"26~S~o~N~N","28~S~o~N~N","30~S~o~N~N","32~S~o~N~N","34~S~o~N~N","36~S~o~N~N"/)   
  res@tmYLMinorValues         = fspan(26,36,26)
  res@tmYLMinorOn             = True  

  res@mpFillOn                = True
  res@mpOutlineOn             = False   
  res@cnFillDrawOrder         = "PreDraw"
  res@mpDataBaseVersion       = "MediumRes"
  res@mpDataSetName           = "Earth..4"
  res@mpAreaMaskingOn         = True
  res@mpMaskAreaSpecifiers    = (/"China","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/)
  res@mpLandFillColor         = "white"
  res@mpInlandWaterFillColor  = "white"
  res@mpOceanFillColor        = "white"
  res@mpOutlineBoundarySets   = "NoBoundaries"

;>--------------------------------------------<
; set for the plot
;>--------------------------------------------<

  res@cnFillOn                = True               
  res@cnLinesOn               = False            
  res@cnLevelSpacingF         = 0.5            
  res@gsnSpreadColors         = True         
  res@lbLabelAutoStride       = True
  res@gsnAddCyclic            = False   
            
  map = gsn_csm_contour_map(wks,final,res)

;>------------------------------------------------------------<
;                      add China map
;>------------------------------------------------------------<

  cnres           = True
  cnres@china     = True        
  cnres@river     = False      
  cnres@province  = True   
  cnres@nanhai    = False      
  cnres@diqu      = True      
  chinamap        = add_china_map(wks,map,cnres)
;>------------------------------------------------------------<
  
  station             = asciiread("/home/lhz/sta_pre_days/station156.txt",(/156,3/),"float")  
  
  res2                 = True                       
  res2@gsMarkerIndex   = 16                       
  res2@gsMarkerSizeF   = 6.               
  res2@gsMarkerColor   = "black"                  
  res2@tfPolyDrawOrder = "PostDraw"  
  res2@cnFillDrawOrder = "PostDraw"  
   
  plots=gsn_add_polymarker(wks,map,station(:,1),station(:,2),res2)
  delete(res2)  
  draw(map)
  frame(wks)
end



数据：   
sta:
56038   98.10   32.98
56079  102.97   33.58
56097  104.25   33.27
56144   98.58   31.80
56146  100.00   31.62
56147   98.83   31.22
56152  100.33   32.28
56158  100.67   31.40
56164  100.98   32.27
56167  101.12   30.98
56168  102.07   31.48
56171  101.70   32.90
56172  102.23   31.90
56173  102.55   32.80
56178  102.35   31.00
56180  103.85   31.68
56181  103.70   30.68
56182  103.60   32.67
56183  103.62   31.50
56184  103.17   31.43
56185  102.98   32.08
56186  104.20   31.33
56187  103.87   30.75
56188  103.67   31.00
56189  103.93   30.98
56190  104.55   31.55
56193  104.52   32.42
56194  104.45   31.63
56195  104.73   31.80
56196  104.73   31.45
56197  104.17   31.15
56198  104.50   31.32
56199  104.68   31.03
56247   99.10   30.00
56251  100.32   30.93
56257  100.27   30.00
56263  101.88   30.88
56267  101.02   30.03
56272  103.88   30.82
56273  102.82   30.38
56276  103.82   30.45
56278  102.77   30.07
56279  102.93   30.15
56280  103.12   30.08
56281  103.50   30.18
56284  103.43   30.45
56285  103.52   30.60
56286  104.25   30.55
56287  103.00   29.98
56288  103.92   30.58
56289  103.87   30.20
56290  104.18   30.78
56291  104.28   30.93
56295  104.55   30.38
56296  104.43   30.85
56297  104.15   30.02
56298  104.60   30.13
56357  100.30   29.05
56371  102.23   29.92
56373  102.85   29.78
56374  101.97   30.05
56376  102.63   29.35
56378  102.35   29.23
56380  103.35   29.88
56381  103.50   30.07
56382  103.60   29.73
56383  103.87   29.83
56384  103.48   29.60
56385  103.33   29.52
56386  103.75   29.57
56387  103.27   29.23
56389  103.95   29.20
56390  104.07   29.67
56391  103.82   30.08
56393  104.85   29.77
56394  104.43   29.45
56395  104.67   29.52
56396  104.77   29.35
56399  104.98   29.18
56441   99.28   28.72
56443   99.80   28.93
56459  101.27   27.93
56462  101.50   29.00
56473  102.77   28.95
56474  102.17   28.55
56475  102.52   28.65
56478  102.43   28.30
56479  102.85   28.00
56480  103.53   28.82
56485  103.58   28.27
56487  103.13   28.33
56490  103.90   28.95
56491  104.57   28.70
56492  104.60   28.80
56493  104.97   28.85
56494  104.15   28.65
56496  105.23   28.32
56498  104.50   28.15
56499  104.78   28.38
56565  101.52   27.43
56569  102.18   27.42
56571  102.27   27.90
56575  102.55   27.37
56578  102.75   27.07
56580  102.80   27.72
56584  103.25   27.70
56592  104.52   28.43
56593  104.92   28.58
56665  101.85   26.68
56666  101.72   26.58
56670  102.12   26.92
56671  102.25   26.65
56674  101.73   26.50
56675  102.58   26.65
57204  105.22   32.57
57206  105.85   32.43
57208  105.52   32.28
57216  106.83   32.35
57217  106.28   32.23
57237  108.03   32.07
57303  105.92   31.73
57304  105.17   31.67
57306  105.97   31.58
57307  105.08   31.10
57308  105.38   31.22
57309  105.88   30.98
57313  106.77   31.87
57314  106.07   31.35
57315  106.40   31.53
57317  106.42   31.03
57318  106.55   31.07
57320  107.22   31.93
57324  107.08   31.58
57326  107.72   31.37
57328  107.50   31.20
57329  107.85   31.10
57401  105.37   30.87
57402  105.70   30.77
57405  105.55   30.50
57407  105.03   30.28
57408  105.35   30.12
57411  106.10   30.78
57413  106.97   30.85
57414  106.45   30.52
57415  106.63   30.52
57416  106.93   30.33
57417  106.28   30.35
57420  107.18   30.75
57503  105.12   29.62
57507  105.30   29.33
57508  105.37   29.15
57600  105.07   28.72
57603  105.83   28.82
57604  105.38   28.78
57605  105.82   28.03
57608  105.43   28.17

baoyu:
0
0.2
0.066666667
0.066666667
0.066666667
0
0
0
0
0
0.066666667
0.066666667
0
0
0
0.066666667
2.2
0
0.066666667
0
0.066666667
3.333333333
2.8
2.266666667
2.4
3.266666667
1.2
3.8
3.333333333
2.2
2.4
2.066666667
1.533333333
0.2
0.066666667
0.066666667
0
0.066666667
2.4
1.2
2.6
3.933333333
2.6
3.866666667
3.666666667
2.8
2.333333333
1.933333333
4.866666667
2.266666667
2.6
2.266666667
1.666666667
2.066666667
1.333333333
2.066666667
2.466666667
0.133333333
0.133333333
2.4
0.133333333
0.733333333
1.066666667
4.133333333
3.533333333
3.933333333
3
3.733333333
4
3.333333333
1
3.133333333
3.066666667
2.533333333
2.533333333
2.866666667
2.266666667
1.866666667
2.666666667
0
0
0.333333333
0
0.6
1.4
0.866666667
1.6
1.466666667
2.466666667
1.2
0.533333333
2.8
3
1.8
2.266666667
2
2.266666667
1.466666667
2.933333333
0.666666667
1.333333333
1.8
1.933333333
1.333333333
0.466666667
1
2.333333333
1.933333333
2.466666667
1.933333333
2.333333333
2.866666667
1.4
2.333333333
2.466666667
2.466666667
3.466666667
3.6
3.933333333
4.066666667
3.333333333
2.666666667
2.8
2.266666667
2.266666667
3
3.266666667
2.333333333
3
2.533333333
3.266666667
4.266666667
3.8
5.066666667
3.4
3.933333333
2.666666667
2.066666667
1.6
2.866666667
3.066666667
2.866666667
3
2.4
2.4
1.8
2.333333333
2.933333333
2.833333333
2.533333333
2.2
2.4
2.333333333
1.866666667
0.866666667
1.866666667

广帅帅

这两天也写了点关于站点插值格点的程序，就楼主这个问题，发表一些粗略的见解。首先明白函数中的rscan( /10,7,4,1/)的中几个参数的意思。NCL 的插值方法采用的是Cressman 和  Barnes，楼主可以抽空看下Cressman插值原理，简单说在插值过程中，会确定插值过程中的影响因子和影响半径，然后先采用大的影响半径进行插值，再用小的影响半径进行继续插值和修订。那么这个rscan( /10,7,4,1/）的几个参数的意思就是前后插值的影响半径。NCL 官网没有对此作太多的说明，只说明括号内最多有四个参数，并且单调递减。根据我对脚本修改的影响来看，这个几个参数就是针对插值的背景场的”格距“，以你的脚本简单来说，你采用的分辨率是0.1*0.1,所以，你现在采用10、7、4、1就是代表第一次插值的影响半径是10除以0.1，也就是100个影响格点。所以，影响范围太大了，因此可能导致数据插值到省外的情况，你可以试试（0.5,0.1,0.05）试试效果会不会好一点~
（以上只是本人理解，如有不对还望见谅，勿责）

淹死的鱼缺氧

多谢楼主分享

just_for_you！

赞。。。。

546243727

请教各位大神：
1.  olon = fspan(97,109,130） ，olat = fspan(26,35,100)  控制了分辨率吗 ， 是把分辨率设置成0.1x0.1了吗？

2.  数据是四川156个站的站点数据，为什么插出省外如此多？ 是 rscan = (/10,7,4,1/)  和
final = obj_anal_ic_deprecated(lon,lat,baoyu,olon,olat,rscan,False)  用得不对吗？

3.  数据没有缺测是不是不用写baoyu@_FillValue = -999.000000？

4. 怎么MARK掉四川以外的地区呢？  

谢谢各位大神！

aas314

学习了   大赞

dannyjoyride

ArcGIS中两三下不就出来了？

timevigor

学习学习，感谢分享

淡如水

  res@mpOceanFillColor = 0     ;用白色填充海洋  0是colormap的索引值

caimint

你好， 请问你解决这个问题了吗？我也遇到了一样的问题，无法mask其他区域谢谢~

546243727

caimint 发表于 2016-1-24 13:23
你好， 请问你解决这个问题了吗？我也遇到了一样的问题，无法mask其他区域谢谢~

;>------------------------------------------------------------<
;            add SiChuan map
;>------------------------------------------------------------<


  filename = "/home/Administrator/SHP/shengjie/bou2_4p.shp"   
   
  f = addfile(filename, "r")   
  
  NAME=(/f->NAME/)  
  asciiwrite ("/home/Administrator/SHP/NAME.txt", NAME)
  sichuan=(/NAME(205)/)  


  geometry = f->geometry  
  segments = f->segments  
  geomDims = dimsizes(geometry)  
  segsDims = dimsizes(segments)  
   

  geom_segIndex = f@geom_segIndex  
  geom_numSegs  = f@geom_numSegs  
  segs_xyzIndex = f@segs_xyzIndex  
  segs_numPnts  = f@segs_numPnts  
  
  lines         = new(segsDims(0),graphic)   
  numFeatures   = geomDims(0)  
  
  
  dlon   = f->x  
  dlat   = f->y  
  
  plres                   = True        
  plres@gsEdgesOn         = True         
  plres@gsEdgeColor       = "black"     
  plres@cnFillDrawOrder   = "PostDraw"  
  plres@tfPolyDrawOrder   = "draw"      
  
  segNum = 0  

  do i=0, numFeatures-1   
         plres@gsFillColor = "gray"
     if( NAME(i).eq.sichuan) then  
         plres@gsFillColor = "transparent"
     end if  
         startSegment = geometry(i, geom_segIndex)      
         numSegments  = geometry(i, geom_numSegs)   
         do seg=startSegment, startSegment+numSegments-1  
            startPT = segments(seg, segs_xyzIndex)      
            endPT = startPT + segments(seg, segs_numPnts)-1
            lines(segNum) = gsn_add_polygon(wks, map, dlon(startPT:endPT), dlat(startPT:endPT), plres)  
            segNum = segNum + 1  
         end do        
  end do  

  delete(plres)  

  draw(map)

  frame(wks)

end