使用NCL绘制安徽省的轮廓图

使用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

绘制的图形如下：