【班门弄斧】用MeteoInfo“制作”各省地图+各省地图下载

beesfly

3q,楼主分享的东西很有用

简玉

初学者飘过，木有金钱下不了各种资源。。。

liujichen

风行前辈新年好！
使用您提供的个别省份shape还有自己制作的shape，会报一些错误：
ERROR 10：Pointer 'hGeom' is NULL in 'OGR_G_GetGeometryCount'.
ERROR 10：Pointer 'hGeom' is NULL in 'OGR_G_GetPointCount'.
……
也不是所有的shape都不行，有的就没问题
很纠结
求您指点！

liujichen

风行前辈新年好！
使用您提供的个别省份shape还有自己制作的shape，会报一些错误：
ERROR 10：Pointer 'hGeom' is NULL in 'OGR_G_GetGeometryCount'.
ERROR 10：Pointer 'hGeom' is NULL in 'OGR_G_GetPointCount'.
……
也不是所有的shape都不行，有的就没问题
很纠结
求您指点！

liujichen

确切地说，好像是在进行某几个省份组合的时候就会出现这样的问题。。。。

MeteoInfo

liujichen 发表于 2016-2-16 10:58
确切地说，好像是在进行某几个省份组合的时候就会出现这样的问题。。。。

说具体些，最好提供出错的文件用于测试

liujichen

MeteoInfo 发表于 2016-2-16 11:08
说具体些，最好提供出错的文件用于测试

最终目的是使用站点数据画南部几个省的年降水分布
使用meteoinfo制作了广东省的shape文件，没有问题，可以画出广东省的年降水分布

1. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   
2. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
   
3. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
   
4. ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
   
5. load "./shapefile_mask_data.ncl"
   
6. 
   
7. begin
   
8. 
   
9. fili = "./nan.txt"
   
10. arg = asciiread(fili,-1,"string")
    
11. count = str_fields_count(arg," ")
    
12. line = dimsizes(arg)
    
13. 
    
14. data = new( (/line,count(0)-1/), string)
    
15. do j = 0, line-1
    
16. do i = 0, count(0)-2
    
17.      data(j,i) = str_get_field(arg(j),i+1," ")
    
18. end do
    
19. end do
    
20. ;print(data)
    
21. station = stringtofloat(data(:,0))
    
22. lat = stringtofloat(data(:,1))
    
23. lon = stringtofloat(data(:,2))
    
24. tem = stringtofloat(data(:,3))
    
25. tem_h = stringtofloat(data(:,4))
    
26. tem_l = stringtofloat(data(:,5))   ;count of rainfall in dali
    
27. pa = stringtofloat(data(:,6))
    
28. rain20 = stringtofloat(data(:,7))
    
29. rn = stringtofloat(data(:,8))*365
    
30. ;wind = stringtofloat(data(:,9))
    
31. 
    
32. ;range of lat(24.6,26.8) and lon(98.8,101.1)
    
33. minlat      =   18.0
    
34. maxlat      =   31.0
    
35. minlon      =   97.0
    
36. maxlon      =  118.0
    
37. olon = new(211,"float")
    
38. olat = new(131,"float")
    
39. data1 = new((/131,211/),"float")
    
40. do i=0,210
    
41. olon(i) = minlon+i*0.1
    
42. end do
    
43. do l=0,130
    
44. olat(l) = minlat+l*0.1
    
45. end do
    
46. 
    
47. olon!0          = "lon"
    
48. olon@long_name  = "lon"
    
49. olon@units      = "degrees-east"
    
50. olon&lon        = olon
    
51. olat!0          = "lat"
    
52. olat@long_name  = "lat"
    
53. olat@units      = "degrees_north"
    
54. olat&lat        = olat
    
55. 
    
56. lon@units       = olon@units
    
57. lat@units       = olat@units
    
58. 
    
59. rn@_FillValue = -9999.0
    
60. ;obj_anal_ic_deprecated--Iterative correction objective analysis (Cressman, Barnes)
    
61. rscan = (/1,.9,.8,.7,.6,.4,.1/)
    
62. data1 = obj_anal_ic_deprecated(lon,lat,rn,olon,olat,rscan, False)
    
63. 
    
64. shp_filename = "./nanbu/test-nan.shp"
    
65. 
    
66. ;county = new(line,"string")
    
67. ;county = (/"Yunlong","Yangbi","Yongping","Dali","Binchuan","Midu",\
    
68. ;          "Xiangyun","Weishan","Jianchuan","Eryuan","Heqing","Nanjian"/)
    
69. 
    
70. data_mask    = shapefile_mask_data(data1,shp_filename,True)
    
71. 
    
72. wks = gsn_open_wks("X11","rain_count_nb")
    
73. gsn_define_colormap(wks,"BlAqGrYeOrRe")
    
74. 
    
75. res                     = True
    
76. res@gsnDraw             = False
    
77. res@gsnFrame            = False
    
78. 
    
79. res@gsnAddCyclic        = False
    
80. ;res@mpDataSetName       = "Earth..4"                                 
    
81. ;res@mpDataBaseVersion   = "MediumRes"
    
82. ;res@mpOutlineOn         = True
    
83. 
    
84. res@mpMinLatF           = minlat
    
85. res@mpMaxLatF           = maxlat
    
86. res@mpMinLonF           = minlon
    
87. res@mpMaxLonF           = maxlon
    
88. 
    
89. res@mpGeophysicalLineThicknessF= 2.
    
90. res@mpNationalLineThicknessF= 2.
    
91. 
    
92. ;res@mpLimitMode = "LatLon"
    
93. ;res@mpLambertParallel1F = .001
    
94. ;res@mpLambertParallel2F = 89.999
    
95. 
    
96. res@cnFillOn      = True
    
97. res@cnLinesOn     = False
    
98. res@cnLineLabelsOn = False
    
99. 
    
100. res@pmTickMarkDisplayMode = "Always"
     
101. 
     
102. map_mask  = gsn_csm_contour_map(wks,data_mask,res)
     
103. lnres                  = True
     
104. lnres@gsLineColor      = "black"
     
105. lnres@gsLineThicknessF = 2.0
     
106. line_mask = gsn_add_shapefile_polylines(wks, map_mask, shp_filename, lnres)
     
107. 
     
108. draw(map_mask)
     
109. 
     
110. ;tres                    =True
     
111. ;tres@txFontHeightF      = 0.01
     
112. ;do i = 0, line-1
     
113. ;    gsn_text(wks,map_mask,county(i),lon(i),lat(i),tres)
     
114. ;end do
     
115. 
     
116. frame(wks)
     
117. 
     
118. end
     

复制代码

。
可是如果使用meteoinfo制作了广东+广西+云南+海南几个省的shape文件，同样的NCL脚背就会报错：
ERROR 10：Pointer 'hGeom' is NULL in 'OGR_G_GetGeometryCount'.
ERROR 10：Pointer 'hGeom' is NULL in 'OGR_G_GetPointCount'.
……
附上脚本

MeteoInfo

liujichen 发表于 2016-2-17 08:59
最终目的是使用站点数据画南部几个省的年降水分布
使用meteoinfo制作了广东省的shape文件，没有问题，可 ...

你是要用shape文件在ncl中做mask功能吗？有可能是ncl的bug（只能用一个polygon图元做mask？），我不是ncl用户也不是很清楚。

siljechen

好东西，即便是现在都好需要的！！

liujichen

MeteoInfo 发表于 2016-2-17 09:26
你是要用shape文件在ncl中做mask功能吗？有可能是ncl的bug（只能用一个polygon图元做mask？），我不是ncl ...

是的，实现mask这么悲剧？