使用NCL绘制安徽省的轮廓图,其它省可以参照设置:
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.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 = 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"/)
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
filename = "/home/PIV/ChinaMap/np/bou2_4l.shp"
;--- Open workstation.
wks = gsn_open_wks("png","Anhui")
;---Create the map
map = create_map(wks,"Anhui of China")
;---Attach the polylines
pres = True
pres@gsLineColor = "blue"
pres@gsLineThicknessF =2
poly = gsn_add_shapefile_polylines(wks,map,filename,pres)
;---Drawing the map will also draw the attached polylines.
draw(map)
frame(wks)
end
图形绘制的结果为:
上面的图形没有添加长江和淮河,下面使用shapefile文件,添加三级以上河流的信息:
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.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 = 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"/)
res2@mpProvincialLineColor="blue"
res2@mpProvincialLineThicknessF =2
;---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
;添加三级以上河流
filename1 = "/home/PIV/ChinaMap/sanjiheliu/hyd2_4l.shp"
filename2 = "/home/PIV/ChinaMap/sanjiheliu/hyd2_4p.shp"
;--- Open workstation.
wks = gsn_open_wks("png","Anhui")
;---Create the map
map = create_map(wks,"Anhui of China")
;---Attach the polylines
pres = True
pres@gsLineColor = "black"
pres@gsLineThicknessF =1
poly1 = gsn_add_shapefile_polylines(wks,map,filename1,pres)
pres2 = True
pres2@gsLineColor = "blue"
pres2@gsLineThicknessF =1
poly2 = gsn_add_shapefile_polygons(wks,map,filename2,pres2)
;---Drawing the map will also draw the attached polylines.
draw(map)
frame(wks)
end
绘制图形如下:
如果想配置更多的图形属性,那么可以这样做:比如只画安徽的轮廓,填充绿色,边界白色,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
绘制图形如下:
下面的代码绘制了安徽气象站点的分布图:
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 = False
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
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")
;*************************************************
; 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
;读入站点位置数据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 = 10. ; polymarker size
res@gsMarkerColor = "blue" ; polymarker color
plots=gsn_add_polymarker(wks,map,sd(:,2),sd(:,1),res);注意经纬度不能错
;
;---Drawing the map will also draw the attached polylines.
draw(map)
frame(wks)
end
绘制的图形如下: