对GMTSAR的D-InSAR得到的结果数据有很多，其中栅格数据类型为grd。
grd数据其实本质上是一种HDF数据（NC数据），在envi或者arcgis是无法打开直接这种的数据。
本文只对los_ll.grd数据进行分析，但分析方法、处理步骤也适用于其他的grd数据。本文分为四个步骤：
1.怎么在python读取grd
2.读取grd的元信息
3.地理变换矩阵
4.grd转换为tif
以上过程是渐进的，实操如下：
怎么在python读取grd
使用gdal可以轻松读取grd，请看下面的代码。
from osgeo import gdal, osrfile = r'los_ll.grd'ds = gdal.Open(file)data = ds.GetRasterBand(1).ReadAsArray()
grd的数据就读取到电脑内存中，而data 是从栅格文件 los_ll.grd 读取的数据数组。它是一个 NumPy 数组，代表了栅格文件中的数值数据。
在pycharm中，data数据如下所示：
data数据的大小为2620x3480,说明data数据的高有2620行，data数据的宽有3480列。
换句话说，2620 表示栅格数据集在垂直方向上有 2620 个像素。这意味着如果你把这个数据集想象成一张图像，那么它的高度是 2620 像素。
同样地，3480 表示栅格数据集在水平方向上的宽度，即从左到右的像素数量。
记住这个2620x3480，后面要用到。
我们可以用matplotlib把这个los_ll.grd文件可视化出来。
import matplotlib.pyplot as pltfrom osgeo import gdal, osrfile = r'los_ll.grd'ds = gdal.Open(file)data = ds.GetRasterBand(1).ReadAsArray()plt.imshow(data),plt.show()
data数据可视化如下所示：
读取grd的元信息
通过gdald的GetMetadata接口，获取grd的元信息。
from osgeo import gdal, osrfile = r'los_ll.grd'ds = gdal.Open(file)info = ds.GetMetadata()
grd的元信息如下所示：
{'lat#actual_range': '{34.9208333333,36.0125}','lat#axis': 'Y','lat#long_name': 'latitude','lat#standard_name': 'latitude','lat#units': 'degrees_north','lon#actual_range': '{-254.005555556,-252.072222222}','lon#axis': 'X','lon#long_name': 'longitude','lon#standard_name': 'longitude','lon#units': 'degrees_east','NC_GLOBAL#Conventions': 'CF-1.7','NC_GLOBAL#description': '','NC_GLOBAL#GMT_version': '6.6.0 [64-bit]','NC_GLOBAL#history': 'gmt xyz2grd llpb -R-254.005555556/-252.072222222/34.9208333333/36.0125 -I2s/1.5s -r -fg -Glos_ll.grd -bi3f','NC_GLOBAL#node_offset': '1','NC_GLOBAL#title': '','z#actual_range': '{-88.73101806640625,51.03330993652344}','z#long_name': 'z','z#_FillValue': 'nan'}
从上图可以知道grd的经纬度的范围。
有意思的是grd的经度范围是-254.005555556,-252.072222222，这明显不符合gdal中的对经度定义的范围。
我们需要对这个小于-180的经度进行改正，即加上360°。让经度范围在-180°到180°之间。
即改正后的经度范围是 105.99444444400001，107.927777778
此外，我们还可以通过经纬度的范围、数据的行列大小，推算出该X 、Y方向上的像素大小(空间分辨率)
X空间分辨率= (107.927777778 -105.99444444400001)/3840 = 0.0005034722223958325
Y空间分辨率= (36.0125-34.9208333333)/2640= 0.0004135101010227273
地理变换矩阵
我们可以构建栅格数据集的地理变换矩阵，这个矩阵提供了栅格数据集与其地理坐标系统之间的空间关系。具体来说，这六个元素按照以下顺序定义：
Transform = (x_origin, x_scale, x_rotation, y_origin,y_rotation, y_scale)六个参数分别为：
左上角x坐标， 水平分辨率，旋转参数， 左上角y坐标，旋转参数，竖直分辨率。一般来说，旋转参数都为0。
x_origin表示栅格左上角像素中心的 X 坐标。y_origin表示栅格左上角像素中心的 Y 坐标。x_scale表示 X 方向上的像素大小（也称为分辨率）。y_scale表示 Y 方向上的像素大小（负数是因为 Y 值随着向下增加而减小）。x_rotation和y_rotation通常为 0，除非栅格被旋转。
现在我们知道了影像的经纬度范围。
我画了一张草图，不用算就可以知道左上角像素中心的 X Y坐标。
x1最小，y1最大。即左上角经度最小、左上角纬度最大
在上面的例子中，los_ll.grd的地理变换矩阵是
（105.99444444400001，0.0005034722223958325，0，36.0125，0，0.0004135101010227273）
grd转换为tif
整合以上内容，一步一步地编写代码，实现grd转换为tif。
暂时把无效值设置为0，全部代码如下：
#!/usr/bin/env python# -*- coding: utf-8 -*-# @File : grd2tiff.py''''''# todofrom osgeo import gdal, osrimport numpy as npimport osimport matplotlib.pyplot as pltclass grd2tiff:def __init__(self, file):self.file = fileself.outfile = os.path.splitext(file)[0] + '.tif'def grd2tiff(self) -> None:"""Do not return anything, modify nums in-place instead."""ds = gdal.Open(self.file)info = ds.GetMetadata()data = ds.GetRasterBand(1).ReadAsArray()# 把nan改成0data[np.isnan(data)] = 0# 84坐标spatialref = osr.SpatialReference()spatialref.ImportFromEPSG(4326)proj = spatialref.ExportToWkt()cols = ds.RasterXSizerows = ds.RasterYSizey_range = [float(info['lat#actual_range'].split(',')[0][1:]), float(info['lat#actual_range'].split(',')[1][:-1])]x_range = [float(info['lon#actual_range'].split(',')[0][1:]), float(info['lon#actual_range'].split(',')[1][:-1])]if x_range[0] < -180:x_range[0] += 360x_range[1] += 360x_res, y_res = (x_range[1] - x_range[0])/cols, (y_range[1] - y_range[0])/rowsleft_x = x_range[0]left_y = y_range[1]Transform = (left_x,x_res,0,left_y,0,-1*y_res)driver = gdal.GetDriverByName('GTiff')self.saveTif( data, cols, rows, driver, proj, Transform)a = 0def saveTif(self, array, cols, rows, driver, proj, Transform):'''@todo single band tif file to save@param array: data array@param cols: cols of data array@param rows: rows of data array@param driver:@param proj: projection@param Transform: coordinate transform@param filename: output filename@return: None'''indexset = driver.Create(self.outfile, cols, rows, 1, gdal.GDT_Float32)indexset.SetGeoTransform(Transform)indexset.SetProjection(proj)Band = indexset.GetRasterBand(1)Band.WriteArray(array, 0, 0)if __name__ == '__main__':file = r'los_ll.grd'a = grd2tiff(file).grd2tiff()
结果可视化如下：
使用qgis转换grd
可以直接使用qgis软件把 grd转换为tif，但是qgis不会把影像的经度范围改正为-180°到180°
所以后续使用矢量去裁剪影像的时候会出错。
把qgis转换的tif，和我们代码转换的tif，放在arcgis可视化如下：
二者相差了360°，所以比例尺很大很大！
自此，完成了对GMTSAR的D-InSAR得到数据进行简单的分析。
如果你还有什么想说的话，可以在留言区留言。
