内容摘要：在上节课的基础上，我们继续学习网格数据的读取和操作。今天我们结合WGM2012的数据，进一步介绍如何采用Geoist的函数来抽取数据，投影变换和画图。
1、读入WGM2012数据
首先，下载WGM2012数据，一般包括：地形、自由空气重力异常、布格重力异常和均衡重力异常四种。默认的是Netcdf的网格文件格式。全球的数据分辨率为2min。上一节课我们讲了GDAL扩展库，如果具体想知道都支持那些格式，可以参考：GDAL支持文件格式列表
在这里我们直接上代码（最重要的是知道GDALGrid这个类）：
import matplotlib.pyplot as plt
from geoist.others.gdal import GDALGrid
from geoist.others.utils import map2DGrid
filename1 = 'D:\\WGM2012_ETOPO1_ponc_2min.grd'
filename2 = 'D:\\WGM2012_Freeair_ponc_2min.grd'
filename3 = 'D:\\WGM2012_Bouguer_ponc_2min.grd'
filename4 = 'D:\\WGM2012_Isostatic_ponc_2min.grd'
gd1, ff = GDALGrid.getFileGeoDict(filename1)
print(gd1)
结果如下：
Bounds: (-180.0000,179.9667,-90.0000,90.0000)
Dims: (0.0333,0.0333)
Shape: (5401,10800)
结果显示这是全球的数据，经度范围-180到179.9667；纬度范围：-90到90。网格间距是0.0333度，网格有5401行，10800列。
这么大的数据其实我没有必要全部都要，所以，可以截取一个区域，方法如下：
gd1.xmin = 70.0
gd1.xmax = 105.0
gd1.ymin = 15.0
gd1.ymax = 45.0
gd1.nx = 800
gd1.ny = 600
# 读取指定区域
topo = GDALGrid.load(filename1, samplegeodict = gd1, resample = True)
fag = GDALGrid.load(filename2, samplegeodict = gd1, resample = True)
bug = GDALGrid.load(filename3, samplegeodict = gd1, resample = True)
iso = GDALGrid.load(filename4, samplegeodict = gd1, resample = True)
注意到上面load函数里面，可以设置samplegeodict的参数，并且要求数据进行重新采样，如果再输入：print(topo._geodict)。
结果如下：
Bounds: (70.0000,105.0000,15.0000,45.0000)
Dims: (0.0438,0.0501)
Shape: (600,800)
可以看到网格间距因为重采样发生了一些变换，但是数据范围和行列都是按照我们设置来的。
那如果，感觉这个区域范围还是太大，再取一块数据怎么办呢？
别着急，可以通过cut函数实现，方法如下：
topoc = topo.cut(80,90, 30,35, align = True)
Cut之后，返回一个新的网格对象topoc，再通过print(topoc._geodict)，结果如下：
Bounds: (80.0313,89.9750,30.0751,34.9833)
Dims: (0.0438,0.0501)
Shape: (99,228)
可以看到由于输入的范围不能正好取到数据，因此，网格边界根据实际情况调整了一下，但是网格间距没变，没有重采样操作。
网格数据读取，差不多就这样子了，下面我们来看看网格数据的投影问题。
2、投影变换
在重磁位场数据处理中，经常要求在直角坐标系下进行计算。因此，需要通过投影变换的方式，将经纬度坐标，变换为直角坐标。
所谓地图投影就是建立地图平面上的点（x,y)和地球表面上的点(x1,y1)之间的函数关系。地图投影中不可避免地存在着变形（长度变形、角度变形），在建立一个投影时不仅要建立（x,y)与(x1,y1)之间的关系，而且要研究投影变形的分布与大小。最常见的投影如图1所示，分为圆柱和圆锥两种。
图1 地图投影样式
首先，我们看一下圆柱投影，比如：图2所示的Mercator圆柱投影，特点是在纬度上面发生变形。
图2 Mercator 圆柱投影
而圆锥投影，如图3的Lambert Conformal Conic projection (LCC)等角圆锥投影， 在经纬度方向都发生一定程度变形。
图3 Lambert 等角圆锥投影
大概了解了投影的概念后，我们看看在python里面怎么做，通过pyproj库，可以使用很多投影，测试一下下面代码：
import pyproj
lcc = pyproj.Proj(proj = 'lcc')
pyproj.Proj.srs
如果安装了pyproj库，应该输出下面信息：
lcc.srs
Out[195]: '+units=m +proj=lcc '
大家注意out后面的字符串，这个就是Proj的string表示形式。也支持EPSG的代码，如：wgs84 = pyproj.Proj(init="epsg:4326")
比如上面两种投影的简写是merc和lcc，大家可能记不住，可以输入：pyproj.pj_list查看全部的定义。
pyproj.pj_list
Out[192]:
{'aea': 'Albers Equal Area',
'aeqd': 'Azimuthal Equidistant',
'airy': 'Airy',
'aitoff': 'Aitoff',
'alsk': 'Mod. Stererographics of Alaska',
'apian': 'Apian Globular I',
'august': 'August Epicycloidal',
'bacon': 'Bacon Globular',
'bipc': 'Bipolar conic of western hemisphere',
'boggs': 'Boggs Eumorphic',
'bonne': 'Bonne (Werner lat_1=90)',
'cass': 'Cassini',
'cc': 'Central Cylindrical',
'cea': 'Equal Area Cylindrical',
'chamb': 'Chamberlin Trimetric',
'collg': 'Collignon',
'comill': 'Compact Miller',
'crast': 'Craster Parabolic (Putnins P4)',
'denoy': 'Denoyer Semi-Elliptical',
'eck1': 'Eckert I',
... ...
实际输出很多，我们就显示到这里。说了这么多，大家要知道正确设置这个string很重要，我们下面分别列出：经纬度、merc和lcc三种的设置方式：
p_jw = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
p_merc = "+proj=merc +lon_0=100 +ellps=WGS84 +datum=WGS84 +no_defs +units=km"
p_lcc = "+proj=lcc +lon_0=100 +lat_0=35 +lat_1=45 +ellps=WGS84 +datum=WGS84 +no_defs +units=km"
这三个string设置好了，就可以进行投影变换了，方法如下：
topo1 = topo.project(p_merc)
topo2 = topo.project(p_lcc)
fig, (ax0, ax1) = plt.subplots(nrows=1, ncols=2,figsize=(12,6))
ax0.imshow(topo1.getData())
ax1.imshow(topo2.getData())
结果如图3所示，原来的地形，变成了两种投影后的样式。
图3 投影后的数据样式，左：Mercator投影，右：LCC投影
如果期望变回来，再转换一次即可，如下形式：
topo0 = topo2.project(p_jw)
我们对比一下投影前后的信息：
print(topo._geodict)
Bounds: (70.0000,105.0000,15.0000,45.0000)
Dims: (0.0438,0.0501)
Shape: (600,800)
print(topo0._geodict)
Bounds: (70.9722,106.3932,13.3594,47.1704)
Dims: (0.0413,0.0413)
Shape: (820,859)
投影转换前后边界会发生一点变化，这时候再重新剪裁一下即可。
3、画几张图
上面介绍完投影功能后，我们再看看怎么把图可以画的好看点，可以通过map2DGrid函数来完成，直接上代码：
from geoist.others.utils import map2DGrid
fig, axes = plt.subplots(nrows=2,ncols=2,figsize=(12,12))
fig.tight_layout()
map2DGrid(axes[0,0],topo,'Topography',10,10, isLeft=True)
map2DGrid(axes[0,1],fag,'Free-air gravity',10,10)
map2DGrid(axes[1,0],bug,'Bouguer gravity',10,10, isLeft=True)
map2DGrid(axes[1,1],iso,'Isostatic gravity',10,10)
结果如图4所示，看看是不是比imshow有点变化了呢？
图4 地形和重力异常网格
一句话结尾：投影变换和网格裁剪都讲完了，有了这些，后面再进行位场计算和反演就可以得心应手了。
网友评论