# -*- coding: utf-8 -*-
# @Time : 2020/10/22 17:16
# @Author : ZhangMin
import os
import pandas as pd
from osgeo import gdal, osr
os.environ['PROJ_LIB'] = r'e:\Anaconda3\envs\cloneTF21\Library\share\proj'
class GRID:
# 读影像文件
def read_img(self, filename):
dataset = gdal.Open(filename) # 打开文件
im_width = dataset.RasterXSize # 栅格矩阵的列数
im_height = dataset.RasterYSize # 栅格矩阵的行数
im_geotrans = dataset.GetGeoTransform() # 仿射矩阵
im_proj = dataset.GetProjection() # 地图投影信息
# 近红外波段
im_data = dataset.ReadAsArray(0, 0, im_width, im_height)
# 将转换为经纬度坐标
# get CRS from dataset
crs = osr.SpatialReference()
crs.ImportFromWkt(dataset.GetProjectionRef())
# create lat/long crs with WGS84 datum
crsGeo = osr.SpatialReference()
crsGeo.ImportFromEPSG(4326) # 4326 is the EPSG id of lat/long crs EPSG 权威编码 OpenGIS WKT坐标系统定义
# get CRS from dataset
t = osr.CoordinateTransformation(crs, crsGeo)
del dataset
return im_data, im_width, im_height, im_geotrans, im_proj, t
if __name__ == '__main__':
run = GRID()
filename = r"L:\2019-2020changjinagkou3\6SWIR_TV\TV.tif"
im_data, im_width, im_height, im_geotrans, im_proj, t = run.read_img(filename)
# (392020.0, 10.0, 0.0, 3499560.0, 0.0, -10.0)
# im_geotrans[1]是像元宽度(影像在宽度上的分辨率);
# im_geotrans[5]是像元高度(影像在高度上的分辨率)
# 这六个参数包括 左上角坐标,像元X、Y方向大小,旋转等信息。 要注意,Y方向的像元大小为负值。
print(im_geotrans)
# 左上角坐标
print(im_geotrans[0])
print(im_geotrans[3])
print(im_proj)
arrSlope = []
lat_long_coor = [] # 用于存储每个像素的(X,Y)坐标
for i in range(im_height):
row = []
for j in range(im_width):
'''
有了下面的公式, 可以把栅格的每个点代入公式求出其在地球上的实际位置
(当然,求出的数值的单位是在坐标系统中定义的那个,而非经纬度坐标, 要求经纬度坐标,还要通过一步转换
'''
px = im_geotrans[0] + i * im_geotrans[1] + j * im_geotrans[2]
py = im_geotrans[3] + i * im_geotrans[4] + j * im_geotrans[5]
# 转换为经纬度坐标
(long, lat, z) = t.TransformPoint(px, py)
# print(long, lat)
lon_lat = [long,lat]
row.append(lon_lat)
lat_long_coor.append(row)
df = pd.DataFrame(lat_long_coor)
df.to_excel("lat_long_coor.xlsx",index=False,header=False)
print(len(lat_long_coor))
python 遥感影像像元坐标读取到excel文件
最新推荐文章于 2023-07-28 15:52:49 发布