内容摘要：网格数据是地球物理数据处理和解释的基础，通常野外测量的离散点也是要先转换为规则的网格数据才能进行进一步的变换与计算。今天我们就谈谈在Geoist中怎么，实现对网格数据的定义和操作接口使用方法。
1、网格数据
（1）基本概念
网格数据这里指的是规则网数据，网格定义如图1所示。
图1 网格文件数据示意图
如果计算机绘制等值线，一般都需要用网格数据来表示。如果不是网格数据，通常需要通过插值算法，先变成网格数据。
在地球物理重磁数据处理和模型计算中，一般也要先插值成规则的网格数据，才能满足FFT变换的数据格式要求。
网格数据定义又很多种，一般都包括文件头和数据体两部分，不同的文件头定义 有差别，因此，常用的包括：ESRI的网格数据定义，Netcdf，Surfer等。
（2）Geoist中的实现
在Geoist的others模块中，实现了对网格数据的定义和操作。基本类的继承关系如下：Dataset->Grid->Grid2D->GMT or GDAL。
而网格信息数据定义类为：GeoDict。
下面我们不急着动手，先看看Geoist中依赖的GDAL库。
2、第三方GDAL库
GDAL(Geospatial Data Abstraction Library)是一个重要的第三方支持库，是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。
有很多著名的GIS类产品都使用了GDAL/OGR库，包括ESRI的ARCGIS，Google Earth和跨平台的GRASS GIS系统。利用GDAL/OGR库，可以直接获得上百种矢量和栅格文件数据的支持。
当然，Python接口也是有的，如果不知道本地已经安装的GDAL版本，可以通过pip show gdal命令查看，结果如下：
(base) D:\Miniconda3\Scripts>pip show gdal
Name: GDAL
Version: 2.3.3
Summary: GDAL: Geospatial Data Abstraction Library
Home-page: http://www.gdal.org
Author: Frank Warmerdam
Author-email: warmerdam@pobox.com
License: MIT
Geoist中，对网格数据类型的格式扩展功能、投影变换、重采样等基于RasterIO库开发，而该库又基于GDAL，所以，这两个第三方库都要支持。
3、一个例子
下面我们结合三个具体的例子，分别说明：
(1)建立一个最简单的网格对象
使用Grid2D和GeoDict两个类，就可以定义一个网格文件，并初始化导入数据。
import matplotlib.pyplot as plt
import numpy as np
from geoist.others.grid2d import Grid2D
from geoist.others.geodict import GeoDict
xmin = 118.5
xmax = 120.5
ymin = 32.0
ymax = 34.0
xdim = 0.25
ydim = 0.25
ncols = len(np.arange(xmin,xmax+xdim,xdim))
nrows = len(np.arange(ymin,ymax+ydim,ydim))
data = np.arange(0,nrows*ncols)
data.shape = (nrows,ncols)
geodict = {'xmin':xmin,
'xmax':xmax,
'ymin':ymin,
'ymax':ymax,
'dx':xdim,
'dy':ydim,
'nx':ncols,
'ny':nrows,}
grid = Grid2D(data,GeoDict(geodict))
plt.figure()
plt.imshow(grid.getData(),interpolation='nearest')
plt.colorbar()
lat,lon = grid.getLatLon(4,4)
print('The coordinates at the center of the grid are %.3f,%.3f' % (lat,lon))
value = grid.getValue(lat,lon)
print('The data value at the center of the grid is %i' % (value))
(2)抽取一个网格内部的区域
有了一个网格对象的示例后，可以使用interpolateToGrid函数对其进行重采样，指定一个新的GeoDict位置后即可，代码如下：
xmin = 119.2
xmax = 119.8
ymin = 32.7
ymax = 33.3
dx = 0.33
dy = 0.33
nx = len(np.arange(xmin,xmax+xdim,xdim))
ny = len(np.arange(ymin,ymax+ydim,ydim))
sdict = {'ymin':32.7,'ymax':33.3,'xmin':119.2,'xmax':119.8,'dx':0.33,'dy':0.33,'nx':nx,'ny':ny}
grid2 = grid.interpolateToGrid(GeoDict(sdict))
plt.figure()
plt.imshow(grid2.getData(),interpolation='nearest'); #'linear', 'cubic', 'nearest'
plt.colorbar();
(3)读取一个实际数据
读取GMT的Netcdf格式网格文件。
from geoist.others.gmt import GMTGrid
gd0, ff = GMTGrid.getFileGeoDict('d:\\out.nc')
grid0 = GMTGrid.load('d:\\out.nc',gd0)
plt.figure()
plt.imshow(grid0.getData())
plt.colorbar()
读取Arcinfo Ascii网格文件。
from geoist.others.gdal import GDALGrid
gd1, ff = GDALGrid.getFileGeoDict('d:\\out-arc.grd')
grid = GDALGrid.load('d:\\out-arc.grd',gd1)
plt.figure()
plt.imshow(grid.getData())
plt.colorbar()
图2 读取后的网格数据
好了，今天就介绍到这里，下一节将继续介绍网格对象的投影变换、内插外插、指定数据提取等扩展功能。
网友评论