分享一个GrADS浙江省县市地图文件

淡如水

附件是浙江边界，里面有很多个浙江边界，有几个是网上下载的，有几个是我根据前面提到的程序做的，你可以分别都作图看一下，记得有两个杭州湾那一块的边界有所区别。其中zhejiang_out是用来白化处理的，但是这种白化处理方法有个问题，就是舟山等岛屿也会被白化掉，所以如果要对浙江进行白化处理，最好用前面的回复提到的另外一种方法，用maskout

淡如水

参考了很多前辈的帖子，同时结合自己的实践经验，前段时间将白化处理的方法进行了总结,与大家分享：
http://blog.sina.com.cn/s/blog_4c72a97e010008ys.html
http://bbs.06climate.com/forum.php?mod=viewthread&tid=13260
主要参考了以上的文章，下面结合自己的实践经验，介绍一下grads白化的处理方法。
主要有两种方法。
1.        首先获得要处理的区域的边界文件（*.bln）,这个文件可以通过以下两种方法获得。在“地图数据取自国家测绘局国家基础地理信息中心1:400万基础地理信息国界数据（http://sms.webmap.cn/default.asp）”网站上获得地区的地图文件（*.shp），将该地图文件用meteoinfo软件打开，点“工具”-“输出地图文件”-“导出surfer  *bln文件”，这样就得到了*.bln文件。如果找不到地区的shp地图文件，也可以在micaps软件的basicGeoInfo文件下找到countryregion.txt文件，该文件包含了全国的行政区域的边界信息，可以用fortran将您需要的地区的信息提取出来（如果是县，直接找到该县，将县的边界信息复制出来另存为xian.bln即可，注意第一行的两个数字分别意味为：经纬度数目、边界外面填白）；有了bln文件以后，一切都好办。接下来用surfer打开bln文件（map-new-basemap的方式打开），然后选择file-export-选择*.mif格式导出，导出的*.mif文件可以用MapInfo打开（表-转入-选择*.mif文件），这时候可以看到地区的地图，如果需要合并地区内的一些县或者区，可以右击一块区域-图层控制-将铅笔符号下面所对应的框打上勾，之后选中您想合并的区域-右键-编辑对象-合并，这样就将区域合并好了，最后将合并后的区域导出，选择表-转出-*.MIF格式，这样就得到了合并区域后的大区域的MIF文件。如果不需要合并地图中的区域，可以直接将surfer导出的*.mif用map_for_grads_region来获得grads地图文件，但是需要根据*.mif的格式对程序做一定的修改（源程序原来是针对mapinfo处理出来的MIF格式的文件的）。接下来打开region_map.for程序，输入MIF的名称，程序处理后就得到了“map_for_grads_region”，将它重名为为*.dat文件，放到grads\dat\文件下面，绘图的时候可以通过'set mpdset *'的方式调用；最后还需要*_out.txt文件通过用province-basemap.gs才能将区域外的信息覆盖。*_out.txt文件需要自己制作，制作该文件必须将要白化的区域的子行政区域合并，得到合并后的bln文件，这个bln文件第一行的第一列数是区域边界点的总数目。以浙江为例，如果要得到这个闭合区域，只能将周边独立的岛删掉才可以。得到的闭合区域命名为zhejiangout.txt,
program main
integer,parameter :: number=3472
real,dimension(number) :: lon,lat
integer ::i
open(unit=1,file='F:\visualfortran\MSDEV98\MyProjects\readlsboundary\zhejiangout.txt')
open(unit=2,file='F:\visualfortran\MSDEV98\MyProjects\readlsboundary\zhejiang_out.txt')
do i=1,number
  read(unit=1,fmt=*)lon(i),lat(i)
end do
write(unit=2,fmt=100)((lon(i),lat(i)),i=1,number)
100 format(3472(f9.5,1x,f8.5,1x))
End
通过上述程序将zhejiangout.txt读成一行输出为zhejiang_out.txt,在该文件开头添加“在该文件前两行添加“117.00000 124.00000  26.00000 33.00000 draw ”，在文件末尾添加“117.00000 26.00000 124.00000 26.00000 124.00000 33.00000 117.00000 33.00000 117.00000 26.00000”，该头分别是绘图区域的经纬度信息，尾部信息要做到头尾相接，也就是开头是点“117.00000 26.00000”，末尾也必须是“117.00000 26.00000”。掌握了这个规律，其实这里的头尾的经纬度信息是可以修改的，但是该经纬度必须是在作图范围内。这样就得到了完整的out文件了。
2.        province-basemap.gs脚本白化的原理是：通过用地图上的点与图形四周四个角的点连线形成闭合区域，然后将每个闭合区域填充为白色；这种方法的致命弱点是：如果某个地区行政区域是由几个不同的闭合区域构成，那么将无法用这种方法实现白化（这是显而易见的）。下面就用第二种方法
3.        打开meteoinfo软件，获得区域的shp文件，分别点击“选择-从左侧框中选择索要导出的图元到右侧框中-输出地图数据-输出数据格式选择grads maskout file-输出，这样就得到了grads maskout文件，该文件由dat文件与ctl（假设该ctl中变量名称为mask）组成。运用该maskout 文件作为背景格点场进行绘图就得到了区域的白化图形。具体调用的时候，首先在gs脚本开头打开地图背景ctl文件（不妨设该ctl是gs脚本中打开的第三个ctl文件），用以下命令：
'define a=oacres(mask.3(t=1),要绘制的气象要素.2,50,25,10)'
'd smth9(maskout(a,mask.3(t=1)))'
这样就得到了白化后的图形。

淡如水

'define a=oacres(mask.3(t=1),要绘制的气象要素.2,50,25,10)'
特别提醒：注意其中插值半径的选择，具体作图的时候可能需要调整插值半径。

淡如水

云知道 发表于 2015-8-2 08:46
怎么没有浙江边界啊

请查看第4页最新发布的回复（该回复的楼上的楼上），是给你的回复，忘了直接点给你回复。

ymcmb

谢谢lz分享

淡如水

有人问起来，就说一说方法

雪颜

淡如水 发表于 2015-8-2 09:53
参考了很多前辈的帖子，同时结合自己的实践经验，前段时间将白化处理的方法进行了总结,与大家分享：
http: ...

如果用maskout的话，插值格点就太少，如果数据间隔大，你把插值半径扩到50，其实误差也是大的

fanrong

刚好用到,太好了

albert007zh

谢谢楼主分享

chenlei1919

谢谢！！！！