使用NCL绘制安徽省的轮廓图 |
您所在的位置：网站首页 › 28个省的轮廓图空白页 › 使用NCL绘制安徽省的轮廓图 |
使用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绘制的图形如下：
|
【本文地址】