Python中从Shapefile提取经纬度及下载LANDSAT影像的方法

1. 确认并转换坐标系 #

Shapefile默认坐标系通常不是经纬度格式（WGS84，EPSG:4326），先检查当前投影，再转换为经纬度坐标系：

import geopandas as gpd

# 读取Shapefile
df_municipalities = gpd.read_file("./CNNNew/Mexico/702825292843_s/mgau2000/agebs_urb_2000.shp")

# 查看当前坐标系
print(df_municipalities.crs)

# 转换为WGS84经纬度坐标系
df_municipalities = df_municipalities.to_crs(epsg=4326)

2. 提取坐标信息 #

提取单个多边形的边界经纬度点 #

转换坐标系后，你之前的代码就能得到标准经纬度点对，调整后更易读：

# 提取每个多边形外边界的（经度, 纬度）列表
df_municipalities['coord_list'] = df_municipalities['geometry'].apply(
    lambda x: [(lon, lat) for lon, lat in x.exterior.coords]
)

获取区域整体边界框（适合影像下载） #

如果目标是下载覆盖所有要素的LANDSAT影像，不需要逐个点，直接提取整体边界范围即可，格式为[最小经度, 最小纬度, 最大经度, 最大纬度]：

# 获取数据集的整体边界框
bbox = df_municipalities.total_bounds
bbox = [bbox[0], bbox[1], bbox[2], bbox[3]]

3. 下载LANDSAT影像 #

用geemap库可以快速对接Google Earth Engine搜索并下载影像，步骤如下：

import geemap

# 创建地图对象
Map = geemap.Map()

# 筛选LANDSAT 8影像（可按需调整年份、时间范围）
landsat_collection = geemap.landsat_timeseries(
    start_year=2020,
    end_year=2020,
    start_date='01-01',
    end_date='12-31',
    roi=bbox
)

# 下载第一景符合条件的影像到本地
geemap.download_ee_image(
    landsat_collection.first(),
    filename='landsat_mexico_region.tif',
    scale=30,  # LANDSAT 8分辨率为30米
    region=bbox
)

注：需先安装geemap（pip install geemap），首次运行需完成Google Earth Engine账号授权。

内容的提问来源于stack exchange，提问作者riskiem