NCL中绘制中国任意省份的精确地图？这个问题我搜索了一下，发现网上都是使用NCL默认的地图做的！
这可了不多，许多中国的领土会不划入争议区！而且仔细对比一下，你会发现NCL给出的行政边界轮廓严重有误！
如果使用国家地理信息网站公布的我国国界和省界的shapefiles文件，完全可以避免上面的问题。
下面的代码是我想突出安徽省的信息，将其标记为白色，可供大家参考：
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"
;----------------------------------------------------------------------
; 绘制安徽省的轮廓图;可以不使用shapefile文件
; 从图中可以看出NCL默认的边界还是很不准的
;----------------------------------------------------------------------
undef("create_map")
function create_map(wks,title)
local a, res2,minlat,maxlat,minlon,maxlon
begin
;---Area to zoom in on.
minlat = 20
maxlat = 40
minlon = 105
maxlon = 125
res2 = True
res2@gsnMaximize = True
res2@gsnDraw = False
res2@gsnFrame = False
res2@mpOutlineOn = True
res2@mpFillOn = False
res2@mpDataBaseVersion = "MediumRes"
res2@mpDataSetName="Earth..4"
res2@mpOutlineSpecifiers=(/"China","Anhui"/) ;NCL自带的地图轮廓，比较粗糙，边界划分失误严重
res2@mpProvincialLineColor="red"
res2@mpProvincialLineThicknessF =4
;---Turn on fancier tickmark labels.
res2@pmTickMarkDisplayMode = "Always"
;---Zoom in on area of interest
res2@mpLimitMode = "LatLon"
res2@mpMinLatF = minlat
res2@mpMaxLatF = maxlat
res2@mpMinLonF = minlon
res2@mpMaxLonF = maxlon
res2@tiMainString = title
;---Create map.
map = gsn_csm_map(wks,res2)
return(map)
end
;---------------------------------------------------------------
begin
;--- Open workstation.
wks = gsn_open_wks("png","Anhui")
;---Create the map
map = create_map(wks,"Anhui of China")
;*************************************************
; Section to add polygons to map.
;*************************************************
filename = "bou2_4p.shp" ;我国公布的国界和省级的Polygon类型的shapefile
f = addfile(filename, "r") ; Open shapefile
NAME=(/f->NAME/)
asciiwrite ("NAME.txt", NAME);从输出的文件中，可以查看第205行显示为"安徽省"，也即NAME(204)
anhui=(/NAME(204)/) ;保存"安徽省"的字符信息,注意strlen(anhui)==6
; anhui=(/"安徽省"/) ;这样定义安徽省，你会发现strlen(anhui)==9
print(f)
print(anhui) ;此处打印"安徽省"的字符会出现乱码，因为NCL不支持宽字符
;
; Read data off shapefile
;
geometry = f->geometry
segments = f->segments
geomDims = dimsizes(geometry)
segsDims = dimsizes(segments)
; Read global attributes
;
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) ; Array to hold polygons
numFeatures = geomDims(0)
lon = f->x
lat = f->y
plres = True ; resources for polylines
plres@gsEdgesOn = True ; draw border around polygons
plres@gsEdgeColor = "white"
colors = (/"blue","green","yellow","red"/)
segNum = 0
do i=0, numFeatures-1
; color assignment (probably a better way to do this?)
if (mod(i,4).eq.0) then
plres@gsFillColor = colors(0)
end if
if (mod(i,4).eq.1) then
plres@gsFillColor = colors(1)
end if
if (mod(i,4).eq.2) then
plres@gsFillColor = colors(2)
end if
if (mod(i,4).eq.3) then
plres@gsFillColor = colors(3)
end if
; 识别是否是安徽的边界
if( NAME(i).eq. anhui) then
plres@gsFillColor = "white" ;安徽省用白色表示
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, lon(startPT:endPT), \
lat(startPT:endPT), plres)
segNum = segNum + 1
end do
end do
;---Drawing the map will also draw the attached polylines.
draw(map)
frame(wks)
end
绘制的结果如下：
代码中使用了一点小的技巧去识别省份的名称，解决了NCL中（中文）宽字符处理的问题。
下面的代码只绘制安徽轮廓地图，绿色填充，黑色边界，并且也用蓝色绘制了NCL自带地图轮廓，从中大家可以看出NCL自带地图局部失真还是较大的：
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"
;----------------------------------------------------------------------
; map_anhui.cnl
; 绘制安徽省的轮廓图;可以不使用shapefile文件
; 从图中可以看出NCL默认的边界还是很不准的
; By Wu Xuping
;----------------------------------------------------------------------
undef("create_map")
function create_map(wks,title)
local a, res2,minlat,maxlat,minlon,maxlon
begin
;---Area to zoom in on.
minlat = 29
maxlat = 35
minlon = 114
maxlon = 120
res2 = True
res2@gsnMaximize = True
res2@gsnDraw = False
res2@gsnFrame = False
res2@mpOutlineOn = True
res2@mpFillOn = False
res2@mpDataBaseVersion = "MediumRes"
res2@mpDataSetName="Earth..4"
res2@mpOutlineSpecifiers=(/"China","Anhui"/) ;NCL自带的地图轮廓，比较粗糙，边界划分失误严重
res2@mpProvincialLineColor="blue"
res2@mpProvincialLineThicknessF =8
;---Turn on fancier tickmark labels.
res2@pmTickMarkDisplayMode = "Always"
;---Zoom in on area of interest
res2@mpLimitMode = "LatLon"
res2@mpMinLatF = minlat
res2@mpMaxLatF = maxlat
res2@mpMinLonF = minlon
res2@mpMaxLonF = maxlon
res2@tiMainString = title
;---Create map.
map = gsn_csm_map(wks,res2)
return(map)
end
;---------------------------------------------------------------
begin
;--- Open workstation.
wks = gsn_open_wks("png","Anhui")
;---Create the map
map = create_map(wks,"Anhui of China")
;*************************************************
; Section to add polygons to map.
;*************************************************
filename = "bou2_4p.shp" ;我国公布的国界和省级的Polygon类型的shapefile
f = addfile(filename, "r") ; Open shapefile
NAME=(/f->NAME/)
asciiwrite ("NAME.txt", NAME);从输出的文件中，可以查看第205行显示为"安徽省"，也即NAME(204)
anhui=(/NAME(204)/) ;保存"安徽省"的字符信息,注意strlen(anhui)==6
; anhui=(/"安徽省"/) ;这样定义安徽省，你会发现strlen(anhui)==9
print(f)
print(anhui) ;此处打印"安徽省"的字符会出现乱码，因为NCL不支持宽字符
;
; Read data off shapefile
;
geometry = f->geometry
segments = f->segments
geomDims = dimsizes(geometry)
segsDims = dimsizes(segments)
; Read global attributes
;
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) ; Array to hold polygons
numFeatures = geomDims(0)
lon = f->x
lat = f->y
plres = True ; resources for polylines
plres@gsEdgesOn = True ; draw border around polygons
plres@gsEdgeColor = "black" ;精确的边界用黑色的表示
segNum = 0
do i=0, numFeatures-1
; 识别是否是安徽的边界
if( NAME(i).eq. anhui) then
plres@gsFillColor = "green" ;安徽省用绿色表示
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, lon(startPT:endPT), \
lat(startPT:endPT), plres)
segNum = segNum + 1
end do
end if
end do
;---Drawing the map will also draw the attached polylines.
draw(map)
frame(wks)
end
图形如下：
屏幕的输出信息如下：
Variable: f
Type: file
filename: bou2_4p
path: bou2_4p.shp
file global attributes:
layer_name : bou2_4p
geometry_type : polygon
geom_segIndex : 0
geom_numSegs : 1
segs_xyzIndex : 0
segs_numPnts : 1
dimensions:
geometry = 2
segments = 2
num_features = 925 // unlimited
num_segments = 978
num_points = 91040
variables:
integer geometry ( num_features, geometry )
integer segments ( num_segments, segments )
double x ( num_points )
double y ( num_points )
double AREA ( num_features )
double PERIMETER ( num_features )
double BOU2_4M_ ( num_features )
double BOU2_4M_ID ( num_features )
integer ADCODE93 ( num_features )
integer ADCODE99 ( num_features )
string NAME ( num_features )
Variable: anhui
Type: string
Total Size: 8 bytes
1 values
Number of Dimensions: 1
Dimensions and sizes: [1]
Coordinates:
(0) 安徽省下面的代码绘制了安徽的气象站点图，此处如果使用NCL默认的地图，你会发现有几个安徽的站点标到了省外去了
。废话少说，看看代码吧：
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"
;----------------------------------------------------------------------
; map_anhui.cnl
; 绘制安徽省的气象站点分布图，并添加一级河流
; By Wu Xuping 2013-11-17
;----------------------------------------------------------------------
undef("create_map")
function create_map(wks,title)
local a, res2,minlat,maxlat,minlon,maxlon
begin
;---Area to zoom in on.
minlat = 29
maxlat = 35
minlon = 114
maxlon = 120
res2 = True
res2@gsnMaximize = True
res2@gsnDraw = False
res2@gsnFrame = False
res2@mpOutlineOn = False
res2@mpFillOn = False
res2@mpDataBaseVersion = "HighRes"
res2@mpDataSetName="Earth..4"
res2@mpOutlineSpecifiers=(/"China","Anhui"/) ;NCL自带的地图轮廓，比较粗糙，边界划分失误严重
res2@mpProvincialLineColor="blue"
res2@mpProvincialLineThicknessF =8
;---Turn on fancier tickmark labels.
; res2@pmTickMarkDisplayMode = "Always"
;---Zoom in on area of interest
res2@mpLimitMode = "LatLon"
res2@mpMinLatF = minlat
res2@mpMaxLatF = maxlat
res2@mpMinLonF = minlon
res2@mpMaxLonF = maxlon
res2@tiMainString = title
res2@tiMainFontHeightF =0.02
res2@tmXBMinorOn = True
;自定义X坐标轴的标签和小标签刻度
res2@tmXBMode = "Explicit"
res2@tmXBValues = (/114,115,116,117,118,119,120/)
res2@tmXBLabels = (/"114~S~o~N~E","115~S~o~N~E","116~S~o~N~E","117~S~o~N~E",\
"118~S~o~N~E","119~S~o~N~E","120~S~o~N~E"/)
res2@tmXBMinorValues = fspan(114,120,31)
res2@tmXBMinorOn = True
;自定义Y坐标轴的标签和小标签刻度
res2@tmYLMode = "Explicit"
res2@tmYLValues = (/29,30,31,32,33,34,35/)
res2@tmYLLabels = (/"29~S~o~N~N","30~S~o~N~N","31~S~o~N~N","32~S~o~N~N",\
"33~S~o~N~N","34~S~o~N~N","35~S~o~N~N"/)
res2@tmYLMinorValues = fspan(29,35,31)
res2@tmYLMinorOn = True
;---Create map.
map = gsn_csm_map(wks,res2)
return(map)
end
;---------------------------------------------------------------
begin
;--- Open workstation.
wks = gsn_open_wks("png","Anhui")
;---Create the map
map = create_map(wks,"Meteorological stations in Anhui")
;*************************************************
;绘制安徽省的地图
;*************************************************
filename = "bou2_4p.shp" ;我国公布的国界和省级的Polygon类型的shapefile
f = addfile(filename, "r") ; Open shapefile
NAME=(/f->NAME/)
asciiwrite ("NAME.txt", NAME);从输出的文件中，可以查看第205行显示为"安徽省"，也即NAME(204)
anhui=(/NAME(204)/) ;保存"安徽省"的字符信息,注意strlen(anhui)==6
; anhui=(/"安徽省"/) ;这样定义安徽省，你会发现strlen(anhui)==9
; print(f)
; print(anhui) ;此处打印"安徽省"的字符会出现乱码，因为NCL不支持宽字符
;
; Read data off shapefile
;
geometry = f->geometry
segments = f->segments
geomDims = dimsizes(geometry)
segsDims = dimsizes(segments)
; Read global attributes
;
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) ; Array to hold polygons
numFeatures = geomDims(0)
lon = f->x
lat = f->y
plres = True ; resources for polylines
plres@gsEdgesOn = True ; draw border around polygons
plres@gsEdgeColor = "black" ;精确的边界用黑色的表示
plres@cnFillDrawOrder = "PostDraw" ; draw polygon first
plres@tfPolyDrawOrder = "draw"
segNum = 0
do i=0, numFeatures-1
plres@gsFillColor = "gray" ;其他省份用灰色表示
; 识别是否是安徽
if( NAME(i).eq. anhui) then
plres@gsFillColor = "green" ;安徽省用绿色表示
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, lon(startPT:endPT), \
lat(startPT:endPT), plres)
segNum = segNum + 1
end do
end do
delete(plres)
;*******************************************************************
;如果有必要可添加一级河流
;*******************************************************************
filename1 = "/d3/SoftWare/ChinaMap/一级河流/hyd1_4l.shp"
filename2 = "/d3/SoftWare/ChinaMap/一级河流/hyd1_4p.shp"
;---Attach the polylines
pres = True
pres@gsLineColor = "red" ;河流的颜色
pres@gsLineThicknessF =2
pres@tfPolyDrawOrder = "PostDraw"
poly1 = gsn_add_shapefile_polylines(wks,map,filename1,pres)
delete(pres)
pres2 = True
pres2@gsEdgesOn = True ; draw border around polygons
pres2@gsEdgeColor = "red" ;河岸的表示
pres2@gsFillColor = "white" ;江湖的颜色
pres2@tfPolyDrawOrder = "PostDraw"
pres2@cnFillDrawOrder = "PostDraw" ; draw polygon first
poly2 = gsn_add_shapefile_polygons(wks,map,filename2,pres2)
delete(pres2)
;*******************************************************************
;读入站点位置数据stationdata:区站号;纬度;经度;观测场海拔(m)
;*******************************************************************
rowsd=81
columnsd=4
sd=asciiread("d.txt",(/rowsd,columnsd/),"float")
;
res = True ; marker mods desired
res@gsMarkerIndex = 16 ; polymarker style
res@gsMarkerSizeF = 15. ; polymarker size
res@gsMarkerColor = "blue" ; polymarker color
res@tfPolyDrawOrder = "PostDraw"
res@cnFillDrawOrder = "PostDraw" ; draw polygon first
plots=gsn_add_polymarker(wks,map,sd(:,2),sd(:,1),res);注意经纬度不能错
delete(res)
;*******************************************************************
draw(map)
frame(wks)
end
绘制出来的图如下：
838
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
