基本上就是一个经纬度转影像坐标的一个操作

之前我用

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 Pointdef 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] #-180yOrigin = geotransform[3]  #90pixelWidth = geotransform[1]pixelHeight = geotransform[5]print(pixelWidth,pixelHeight)from pyproj import Proj# 首先定义要转换的投影坐标系proj1 = Proj("epsg:3413")  # coord_x = lon# coord_y = latcoord_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.RasterXSizenum_rows = raster.RasterYSizeif loc_x<0 or loc_y<0 or loc_x>=num_columns or loc_y>=num_rows:return '999',999data_value = raster.GetRasterBand(1).ReadAsArray(loc_x, loc_y, 1, 1)[0, 0]return data_value