【python + GIS】提取遥感影像经纬度为二维数组

一、应用场景

最近在求非平坦地区的太阳高度角，其中有一个参数就是纬度，因此就想到了把纬度作为值来创建二维数组的想法。

二、参考文献

《读取遥感影像的信息》
《使用python读取tiff文件中的经纬度，并将数据以excel表的形式输出（详细步骤）》

三、代码展示

import os
import numpy as np
from osgeo import gdal
import glob
import time
import math
import copy


DEM_Path = "D:\\DATA\\SEBAL\\0. InitialData\\DEM1984.tif"  #DEM影像
dataset = gdal.Open(DEM_Path)  # 打开tif
# 获取行数列数和地理信息
# geo_information(0):左上像素左上角的x坐标。    
# geo_information(1):w - e像素分辨率 / 像素宽度。    
# geo_information(2):行旋转（通常为零）。
# geo_information(3):左上像素左上角的y坐标。
# geo_information(4):列旋转（通常为零）。
# geo_information(5):n - s像素分辨率 / 像素高度（北半球上图像为负值）   
width = dataset.RasterXSize  # 宽度
height = dataset.RasterYSize  # 高度
geo_information = dataset.GetGeoTransform()
proj = dataset.GetProjection() #地图投影
filevalue= np.array(dataset.ReadAsArray(0, 0, width, height), dtype=float)  # 将数据写成数组，对应栅格矩阵
LongitudeValue=np.zeros((height,width))  #创建0矩阵
LatitudeValue=np.zeros((height,width))
for x in range(height):  # 行     
   for y in range(width):  # 列
    # 赋值经纬度
       LongitudeValue[x][y]=geo_information[0] + y * geo_information[1] + x * geo_information[2]
       LatitudeValue[x][y]=geo_information[3] + y * geo_information[4] + x * geo_information[5]

最后得到的数组和图像分别是这样的：
经度：

纬度：

四、反思与注意事项

1、遥感影像的长度和宽度以及二维数组的行列要对应正确，比如for循环的时候，如果写反了会出现长度不够的错误
2、 创建0矩阵的时候也是要注意参数中的行列对应，如果不行就交换一下
3、这个算法的核心是掌握地理信息的六个参数：
geo_information(0):左上像素左上角的x坐标。
geo_information(1):w - e像素分辨率 / 像素宽度。
geo_information(2):行旋转（通常为零）。
geo_information(3):左上像素左上角的y坐标。
geo_information(4):列旋转（通常为零）。
geo_information(5):n - s像素分辨率 / 像素高度（北半球上图像为负值）