39 lines
1.3 KiB
Python
39 lines
1.3 KiB
Python
import rasterio
|
|
from rasterio.warp import transform
|
|
import matplotlib.pyplot as plt
|
|
import glob
|
|
import os
|
|
|
|
def get_bounds_lonlat(dataset):
|
|
# 直接用元数据获取边界
|
|
left, bottom, right, top = dataset.bounds
|
|
# 四角顺序：左下、左上、右上、右下
|
|
xs = [left, left, right, right]
|
|
ys = [bottom, top, top, bottom]
|
|
# 转为经纬度
|
|
lons, lats = transform(dataset.crs, 'EPSG:4326', xs, ys)
|
|
return list(zip(lons, lats))
|
|
|
|
if __name__ == "__main__":
|
|
tif_dir = "test_geosat1_imgs"
|
|
pattern = os.path.join(tif_dir, "Jilin*.tif")
|
|
all_bounds = []
|
|
for tif in glob.glob(pattern):
|
|
with rasterio.open(tif) as ds:
|
|
corners_lonlat = get_bounds_lonlat(ds)
|
|
all_bounds.append((os.path.basename(tif), corners_lonlat))
|
|
|
|
plt.figure(figsize=(10, 10))
|
|
for name, corners in all_bounds:
|
|
# 闭合多边形
|
|
corners_closed = corners + [corners[0]]
|
|
lons, lats = zip(*corners_closed)
|
|
plt.plot(lons, lats, marker='o', label=name)
|
|
for lon, lat in corners:
|
|
plt.text(lon, lat, f'({lon:.4f},{lat:.4f})', color='red', fontsize=8)
|
|
plt.xlabel("Longitude")
|
|
plt.ylabel("Latitude")
|
|
plt.title("Jilin* 影像边界经纬度分布")
|
|
plt.legend()
|
|
plt.grid(True)
|
|
plt.show() |