基本上就是一个经纬度转影像坐标的一个操作
之前我用
xOrigin = geotransform[0] #-180
yOrigin = geotransform[3] #90
这两个读取出来的分别就是经度和纬度,但是读取极投影为3413的影像时,读取出来的时投影坐标,因此在程序中多做了一步变换
from osgeo import gdal
import os
import pandas as pd
from datetime import datetime, timedelta
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
import numpy as np
# import rasterio
from rasterio.transform import from_origin
from shapely.geometry import Point
def get_location_data(lon,lat,tif_files):
size=625 #12.5km中心半径
f=tif_files
# print(f)
tif_name = os.path.basename(f).split('_')[0]
raster: gdal.Dataset = gdal.Open(f)
geotransform = raster.GetGeoTransform()
#获取栅格影像的左上角起始坐标,像元大小
xOrigin = geotransform[0] #-180
yOrigin = geotransform[3] #90
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
print(pixelWidth,pixelHeight)
from pyproj import Proj
# 首先定义要转换的投影坐标系
proj1 = Proj("epsg:3413")
# coord_x = lon
# coord_y = lat
coord_x,coord_y = proj1(lon,lat)
#主要思路就是计算该坐标与该tif起始坐标差了多少行和列
loc_x = int((float(coord_x) - xOrigin) / pixelWidth)
loc_y = int((float(coord_y) - yOrigin) / pixelHeight)
# print(loc_y)
#知道了多少行和列,就直接读这个行列对应数像元的数值大小,并把读到的数值追加到data这个空数组里面
# with rasterio.open(f) as src:
# # 获取图像的地理转换信息
# transform = src.transform
# # 创建经纬度点的Shapely几何对象
# point = Point(lon, lat)
# # 将经纬度点转换为图像坐标
# lon, lat = point.x, point.y # 重新赋值,确保点在图像范围内
# loc_x, loc_y = ~transform * (lon, lat) # 逆变换,获取图像坐标
num_columns = raster.RasterXSize
num_rows = raster.RasterYSize
if loc_x<0 or loc_y<0 or loc_x>=num_columns or loc_y>=num_rows:
return '999',999
data_value = raster.GetRasterBand(1).ReadAsArray(loc_x, loc_y, 1, 1)[0, 0]
return data_value