给定经纬度坐标或者投影坐标，提取对应位置上栅格的值(python)

在给定一组栅格数据以及给定很多坐标点时，需要提取对应坐标点上的栅格值。如下图所示

 

arcgis中有多值提取到点的操作，也可以使用arcpy进行提取。但是这里存在两个问题。（1）提取完成后需要从arcmap属性表导出为excel文件；（2）如果栅格数据过多，站点数据又很多的时候，这种时候很麻烦，也很崩溃。所以我们打算用python进行处理。

关于arcpy的可以参考Python地理数据处理 十九：arcpy批量处理数据之多值提取至点_extractmultivaluestopoints 批处理_Amyniez的博客-优快云博客

 今天我们使用python中的rasterio和geopandas两个库进行快速提取，最后按照站点名导出所需数据。

注意：我们需要将经纬度数据转为shp文件，这可以在arcmap中实现。

具体代码：

# coding = utf-8

'''

author: max
date: 2023.7.5

# description
这个程序用于提取给定坐标点下（可以是多个点坐标，需要使用shp格式）对应的栅格数值


'''

import os
import rasterio
import geopandas as gpd
import pandas as pd
# import matplotlib.pyplot as plt
from  glob import glob
# from rasterio.plot import show


# 输入路径
input_shp_point = r'C:\Users\15845\Desktop\antarctica_validation_point_projection.shp'
input_raster_path = r'H:\third_paper\antartica\day'
output_extract_path = os.path.join(input_raster_path,'result')
if not os.path.exists(output_extract_path):
    os.makedirs(output_extract_path)


# (1)打开shp文件路径
point_dataset = gpd.read_file(input_shp_point)
point_name = [] # 定义一个列表存储最后的站点名称
# print(point_dataset)
# print(type(point_dataset.Name))
for name in point_dataset.Name:
    point_name.append(name)
print(point_name)



# (2) 从point_dataset提取 x,y,存储到列表中
lon_point = []
lat_point = []
for point in point_dataset['geometry']:
    # print(point.name)
    # print(point.xy[0][0],point.xy[1][0])
    # print(point)
    lon = point.xy[0][0]
    lat = point.xy[1][0]
    lon_point.append(lon)
    lat_point.append(lat)

# 这里站点的名称数和站点个数必须要一致对应上
print(len(point_name))
print(len(lon_point))


# (3)打开待提取的栅格数据路径
raster_filelist = glob(os.path.join(input_raster_path,'*.tif'))
# print(raster_filelist)


# # (4)定义一个字典来存储最后的值
# lst_value = []
# lst_value_total = []
# lst_name_total = []
# for name in point_name:
#     lst_name_total.append(name)
#     lst_value_total.append(lst_value)

# dict_name_lst = dict(zip(lst_name_total,lst_value_total))
# print(dict_name_lst)

# 定义一个嵌套列表
list_lst_value_total = []
for i in range(len(point_name)):
    list_value = []
    list_lst_value_total.append(list_value)
print(list_lst_value_total)

year_month_list = []
# (5)遍历每一个栅格
for raster_file in raster_filelist:
    # print(raster_file)

    # 获取栅格的日期
    year_month = os.path.basename(raster_file).split('.')[0]
    year_month_list.append(year_month)
    raster_file_open = rasterio.open(raster_file)
    print(raster_file_open)
    print((raster_file_open.read(1)))
    # print(raster_file_open.crs)
    # print(raster_file_open.count)

    # 对每个站点进行处理
    if len(point_name) == len(lon_point):
         for i in range(len(point_name)):
             
             # 对每个站点进行处理,获取每个站点的名称，站点经纬度
            point_name_single = point_name[i]
            print('point_name_single',point_name_single)
            point_lon_single = lon_point[i]
            point_lat_single = lat_point[i]
            # print(point_lon_single,point_lat_single)
            #获取对应站点的行列号
            row,col = raster_file_open.index(point_lon_single,point_lat_single)
            print(row,col)
            lst_value = raster_file_open.read(1)[row,col]
            lst_value_new = round(lst_value,2) # 保留两位小数
            print(lst_value_new)

            list_lst_value_total[i].append(lst_value_new)

print(year_month_list)
print(list_lst_value_total[24])


# 将最终的结果输出为csv文件，按照给定的坐标名称

for i in range(len(point_name)):
    point_name_value = point_name[i] # 站点名称
    output_extract_path_file = os.path.join(output_extract_path, "{}.csv".format(point_name_value)) # 设置输出路径和文件名
    
    output_extract_value = list_lst_value_total[i]

    # output,设置为dataFrame
    output_extract_pd = {
        'year_month':year_month_list,
        'lst':output_extract_value
    }

    pd_output_file = pd.DataFrame(output_extract_pd)
    pd.DataFrame.to_csv(pd_output_file, output_extract_path_file, index=False)