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

2024-06-27 00:55| 来源: 网络整理| 查看: 265

NCL中绘制中国任意省份的精确地图？这个问题我搜索了一下，发现网上都是使用NCL默认的地图做的！

这可了不多，许多中国的领土会不划入争议区！而且仔细对比一下，你会发现NCL给出的行政边界轮廓严重有误！

如果使用国家地理信息网站公布的我国国界和省界的shapefiles文件，完全可以避免上面的问题。

下面的代码是我想突出安徽省的信息，将其标记为白色，可供大家参考：

代码中使用了一点小的技巧去识别省份的名称，解决了NCL中（中文）宽字符处理的问题。

下面的代码只绘制安徽轮廓地图，绿色填充，黑色边界，并且也用蓝色绘制了NCL自带地图轮廓，从中大家可以看出NCL自带地图局部失真还是较大的：

屏幕的输出信息如下：

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 绘制出来的图如下：

【本文地址】

公司简介

联系我们