本文主要是介绍python 已知经纬度获取影像DN值或位置，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
基本上就是一个经纬度转影像坐标的一个操作
之前我用
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这篇关于python 已知经纬度获取影像DN值或位置的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
