python画出WPS中namelist的范围

一水天下

最近需要处理WRF Hydro的DEM数据，因此，需要知道流域的精确经纬度范围。WRF Hydro的耦合模式计算了非常大，为了节约计算资源，希望内层嵌套尽量逼近流域边界。但是好像现有的NCL没有这么详细的展示WPS Domain的脚本。顺便说一句，NCL用起来实在是太慢了，虽然气象绘图功能的确是最强的（本人横向对比过R、matlab、Grads）。Python在加入了PyNIO和PyNGL库以后，NCL的基础功能都有了，同时还有很强大的计算库、统计库之类，WingIDE的交互界面也还可以，debug起来也还算看得清楚，虽然交互界面和R Studio和matlab比起来是渣渣，但是比NCL好多了。
所以开始转战python，python也的确没让我失望。
在介绍我的库之前，先说明下。PyNGL库的开发好像要停止了，和开发者联系了，说他们的项目经费被停了。。。。OTZ。不过后面有wrf_python的函数库正在开发，据说比PyNGL要更加稳定。
在使用脚本前。先要准备下面的函数库：
Nio －NETCDF
Ngl －NCL
（特别注意！！！！！！！！，如果你是在2016年11月以前下的pyNGL库，那是不行的，因为里面有bug，11月以后NCAR修复了。就可以用了。）
numpy －数学运算
f90nml －处理f90 namelist

shapefile －处理GIS shapefile地形文件
site －内置

Python版本2x，暂时使用的2.7，2.6没试过。3x肯定不行。2.4以前好像也肯定不行。

操作系统mac OSX，windows不知道。

namelist要自己写。需要注意一点，namelist要写的标准一点，因为程序主要是通过读取namelist中的变量来画图的，如果namelist很糟糕，那就没用了。

精细化的中国省、县shapefile，在http://www.gadm.org/country/上面下载，有其他国家的。中国的部分在附件中有上传。

效果图如下：

                               
登录/注册后可看大图

                               
登录/注册后可看大图

主要的特点是在每个domain的四个角上都有很精确的经纬度坐标。实际应用中，在ArcGIS抠出流域以后，使用WGS84坐标看流域边界的最大最小经纬度，通过调namelist，run我的脚本，画出domain。查看内部domian是不是在流域内。（需要注意一点，WRF虽然使用的圆球坐标系，但是它使用的地形资料都是WGS84坐标下的USGS等数据。并且WRF自己没有做任何坐标转换。这样做其实是可以的，具体解释起来比较麻烦。总之原因是，当所有东西都建立在错错误的框架下的时候，那就错错错错的对了。当然，在ArcGIS中使用WGS84以外的坐标的数据是不行的，要转换成WGS84才能使用我的脚本和放入WRF。）

图中可以改动的内容包括shapefile的精度、shapefile的线条特征、坐标投影（目前只有WRF中1、2、3种）、"D01"字体颜色等、还有地图地图底图的填充色等。这些的修改需要对PyNGL和NCL有一定的了解。

使用步骤：
#找python site-packages 库:
import site; site.getsitepackages()
#将库文件wps_dom.py放入site-packages目录下。
#在site-packages下新建文件夹wps_maps, 把shapefile文件都放入wps_maps下。
#写namelist.wps
#修改脚本setting.py中namelist路径和出图路径，当然画图属性也可以在这里设置。


如果你要修改D01的坐标边框的等属性，修改wrf_wps_dom中的opts_arg.如果你要修改D02和D0n的边框颜色粗细和标注的等属性，就要注意了。你要使用到双层嵌套的Ngl.Resource()说明。如下：


tx_res_1                = Ngl.Resources()
#tx_res_1中包含每个domain的字体resource设置。
#分别封装在tx_res_1.D01，tx_res_1.D02.。。。中。
#D01是空的，因为opts_arg可以设置。
#同时每个tx_res_1.D0n都是一个resource说明。需要单独设置。
#ln_res_1相同

tx_res_1.D01                = Ngl.Resources()
tx_res_1.D02                = Ngl.Resources()
tx_res_1.D02.txJust        = "BottomLeft"
ln_res_1                = Ngl.Resources()
ln_res_1.D01                = Ngl.Resources()
ln_res_1.D02                = Ngl.Resources()
ln_res_1.D02.gsLineColor="green"



#tx_res代表标注属性。
#ln_res标示每个domain的外框属性。

修改shapefile也不一样：
#shp_res代表shapefile的属性。
shp_res=True代表使用默认
shp_res=False代表不叠加地图（ln_res和tx_res没有此选项）


shp_res=Ngl.Resource()
shp_.......so on代表使用自己设置属性。
#shp_res和因为是整个的地图，所以只有一层resource。
#shp_res的地图式样有3种精度，通过shp_res.name来设置，3个分别为：CHN_adm1，CHN_adm2，CHN_adm3。3精度最高。例子如下：
shp_res =Ngl.Resource()
shp_res.name="CHN_adm1"





更多其它内容加入和bug修复后面会持续更新。