版權聲明：本文為部落客原創文章，遵循 CC 4.0 BY-SA 版權協定，轉載請附上原文出處連結和本聲明。
平時在做遙感應用的時候經常會拿監測站點資料和遙感影像進行對比，然後做一些拟合什麼的，本部落格就是實作從有經緯度資料的excel/csv表格中提取相應影像的像元值。
示例表格資料：
| Id | X | Y |
| 1 | 108.98054534900 | 34.49123985970 |
| 2 | 109.51830353600 | 36.41146489150 |
| 3 | 110.01799034700 | 33.55611168450 |
| 4 | 107.88837274700 | 32.75899224750 |
| 5 | 107.78129700200 | 33.97251736050 |
| 6 | 109.73245502700 | 37.55360617420 |
| 7 | 108.25718920300 | 37.17289241330 |
| 8 | 109.09000055500 | 33.68698203980 |
| 9 | 109.92320848400 | 35.62469610990 |
| 10 | 107.93338421800 | 34.81270504160 |
| 11 | 110.04813018700 | 35.53844064840 |
| 12 | 108.30914076800 | 33.82205610950 |
| 13 | 108.87783194900 | 34.01241299000 |
Python實作過程：
# encoding: utf-8
# author：三十二号星期八
# date ：2021/2/1
import gdal
import os
import pandas as pd
def get_location_data(location,tif_files):
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]
yOrigin = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
#建立一個空的數組，把for循環裡面得到的值放到這個數組裡面
data = []
for index, row in location.iterrows():
long_value = row["X"]
lat_value = row["Y"]
print(lat_value)
print(yOrigin)
coord_x = long_value
coord_y = lat_value
#主要思路就是計算該坐标與該tif起始坐标差了多少行和列
loc_x = int((float(coord_x) - xOrigin) / pixelWidth)
loc_y = int((float(coord_y) - yOrigin) / pixelHeight)
print(loc_y)
#知道了多少行和列，就直接讀這個行列對應數像元的數值大小，并把讀到的數值追加到data這個空數組裡面
data_value = raster.GetRasterBand(1).ReadAsArray(loc_x, loc_y, 1, 1)[0, 0]
data.append(data_value)
location['NDVI_'+tif_name] = data
print(location)
return location
if __name__ == '__main__':
excel_data = pd.read_excel(r'G:\0test\test\xj.xlsx')
test_tif =r'G:\0test\test\xj.tif'
out_data = get_location_data(excel_data,test_tif)
#輸出為另一個表格
xls_name = r'G:\0test\test\location_value.xlsx'
out_data.to_excel(xls_name)
IDL實作過程：
ps： IDL正常讀取excel資料較難，是以我把測試的excel轉成了csv格式，便于讀取。
;encoding:GB2312
;author:YX
;date:2021-2-1
pro read_csv_value
compile_opt idl2
e=envi(/headless)
fn='G:\0test\test\xj.csv'
raster=e.openraster('G:\0test\test\xj.tif')
raster_data=raster.GetData()
data = read_csv(fn,count=nl,header=header)
LeftPoint_arr=raster.SPATIALREF.TIE_POINT_MAP
PIXEL_SIZE=raster.SPATIALREF.PIXEL_SIZE
minX=LeftPoint_arr[0]
maxY=LeftPoint_arr[1]
;同樣建立一個空數組，存放資料
anglevalue_arr=make_array(4,nl)
for i=0,nl-1 do begin
;讀取csv裡面經緯度數值
loc_X=data.field2[i]
loc_Y=data.field3[i]
;計算該坐标對應的行列号
sample=floor((loc_X-minX)/PIXEL_SIZE[0])
line=floor((maxY-loc_Y)/PIXEL_SIZE[1])
;擷取該行列号對應的栅格值
anglevalue=raster_data[sample,line]
;把經緯度和對應的栅格值放到上面建立的三維數組裡面
anglevalue_arr[0,i]=i
anglevalue_arr[1,i]=loc_X
anglevalue_arr[2,i]=loc_Y
anglevalue_arr[3,i]=anglevalue
endfor
;定義輸出csv的header
header=['ID','X','Y','xj_tif_value']
;輸出
write_csv,'G:\0test\test\IDL_location.csv',anglevalue_arr,header=header
end
結果，左邊為python結果，右邊為IDL結果
在這個過程中發現了兩個小細節：
第一個小細節是計算該經緯度對應那一列的時候，在python裡面是
loc_y = int((float(coord_y) - yOrigin) / pixelHeight)
在IDL裡面為：
line=floor((maxY-loc_Y)/PIXEL_SIZE[1])
python裡面讀取列方向的像元大小時，讀出來是個負值，是以必須用低緯減去高緯（影像左上角的起始緯度），而在IDL裡面讀取列方向的像元大小時，讀出來是個正值，是以要影像左上角的起始緯度減去表格的緯度，也就是高緯減低緯。
第二個小細節，就如下圖所示，arcmap_value那一列字段的值是在arcmap裡面用Spatial analyst——提取分析——值提取至點，這個工具提取出來的，但提取出來的值和IDL和Python裡面的有幾個不一樣，右圖是在arcmap裡面識别出來的值，又和代碼提出來的是一樣樣的。。。
是以，有沒有人能解答一下這個疑惑呢？抱拳了~