NCL中绘制中国任意省份的精确地图

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

绘制出来的图如下：