IDL&python学习——实现根据有经纬度坐标的excel/csv表格提取相应影像的像元值

版权声明：本文为博主原创文章，遵循 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里面识别出来的值，又和代码提出来的是一样样的。。。

                                          

所以，有没有人能解答一下这个疑惑呢？抱拳了~