内容摘要:上次写完Geoist的Flex的岩石圈有效弹性厚度估计功能之后,大家反映描述的例子不好,分明是用现成的例子来敷衍交差。为了表明动手能力,小哥周末从头开始搭建了一个Linux上的python3环境,亲自来趟坑。下面就来看看这次实战的成果。
1、搭建一个虚拟机
小哥早晨起来,在服务器上找了一个虚拟的Ubuntu系统,远程ssh上去,配置了一个python3.7+jupyter Lab的组合。
这样在本地就可以直接用浏览器开始工作了,配置一个能运行geoist的环境,需要安装不少依赖包,这里具体过程就略去了,安装过程中稍微复杂的,是GDAL库,如果编译提示gdal-config类错误,记得安装libgdal-dev,如下命令:sudo apt-get install gdal-bin libgdal-dev
另外,就是为了画地图,需要basemap支持和离线地图数据包。通过,升级matplotlib解决。不行的就试试sudo apt-get install python-mpltoolkits.basemap。
其它的没遇到太多问题,pip3 install ***,一般都可以解决。
安装完jupyter lab后,需要修改配置文件,给个密码,并且设置外网可访问,即可。整个过程大概花了1个多小时,配置结果如图1所示。
ps:内网的用户直接可以访问上面地址,需要密码可以找小哥要
2、数据配置
WGM2012我们已经提到好几回了,全球数据其实就四个grd文件。每个200多M,分别为地形,自由空气重力异常,布格重力异常,还有均衡异常。我们都传到DATA目录下,然后,想研究哪里,就切哪里即可。图1上面的代码我们切出来一块,地形布格异常如图2所示。
原来数据是经纬度坐标,图2中的数据做了一个投影,变为直角坐标,显示如下:
看到图中的黑线了吗?这是海陆边界,大陆上还有一些湖泊,这些线可以作为参考,用以检查由于投影进行仿射变换前后,数据重采样后是否发生畸变。
具体计算过程,代码就不说了,下面是30个格点作为一个区域,开窗计算后的和的误差。
有了计算面上的计算结果,我们也可以区任何一个区域,看看观测结果与计算模型之间的差异,如图6所示,其中橙色曲线为最优估计模型的拟合结果,蓝色为实际数据观测结果。
在模型中有两个不确定性参数,上节课也提到了和,这里面我们用贝叶斯方法,引入pyMC3的MCMC算法来估计这两个模型参数,如图7所示。
需要再啰嗦一句的是,pymc3实现的mcmc算法,在linux上的进程通信,好于win。有些莫名其妙的问题在linux一下子过去了。程序调好,部署到服务器上,还是比较省事。但是没有图形界面,调试效率低了些,最好是都掌握了,相互调用吧。
一句话总结:综上所述,小哥全部测试了一遍Flex模块的估计功能,因为是从开源项目中修改过来的,实践中往往会遇到一些问题,确实小哥也发现了,有的与依赖库版本有关,有的与操作系统有关,具体我就不唠叨了,还是那个原则,有坑我来趟,大家直接走即可。好了,今天就说到这,上面说的都可用,你要感兴趣,遇到具体问题,欢迎找我讨论。