对GMTSAR的D-InSAR得到数据进行简单的分析

ytkz2024-08-202024-09-13

对GMTSAR的D-InSAR得到的结果数据有很多，其中栅格数据类型为grd。

grd数据其实本质上是一种HDF数据（NC数据），在envi或者arcgis是无法打开直接这种的数据。

怎么在python读取grd

使用gdal可以轻松读取grd，请看下面的代码。

from osgeo import gdal, osr
file = 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 plt
from osgeo import gdal, osr
file = 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, osr
file = 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_rotation, y_origin, y_scale)

其中：

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 
'''

'''
# todo
from osgeo import gdal, osr
import numpy as np
import os
import matplotlib.pyplot as plt
class grd2tiff:
    def __init__(self, file):
        self.file = file
        self.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改成0
        data[np.isnan(data)] = 0

        # 84坐标
        spatialref = osr.SpatialReference()
        spatialref.ImportFromEPSG(4326)
        proj = spatialref.ExportToWkt()
        cols = ds.RasterXSize
        rows = ds.RasterYSize

        y_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] += 360
            x_range[1] += 360
        x_res, y_res = (x_range[1] - x_range[0])/cols, (y_range[1] - y_range[0])/rows

        left_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 = 0
    def 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()