使用python以DEM为输入栅格数据，按1000 m的间隔提取海拔等值线

在Python中，我们可以使用一些地理信息处理库，如Rasterio、GeoPandas和matplotlib等来提取DEM（数字高程模型）的海拔等值线。以下是一个简单的步骤概述：

1. **安装所需库**:
首先，确保已经安装了`rasterio`, `geopandas`, 和 `matplotlib`。如果没有，可以分别通过pip安装：

   pip install rasterio geopandas matplotlib

2. **导入模块**:

   import rasterio
   from rasterio.plot import show
   import geopandas as gpd
   import numpy as np

3. **读取DEM数据**:
使用`rasterio`打开DEM文件，并读取其内容：

   with rasterio.open('dem_file.tif') as src:
       dem_data = src.read(1)  # 第1带通常包含海拔信息
       transform = src.transform

4. **计算等值线**:
创建等值线列表，例如每1000米一个等高线：

   elevations = np.arange(src.nodata, src nodata + (src.height * src.width), 1000)

5. **提取等值线**:
使用`scipy`库中的`griddata`函数将DEM数据插值到等值线位置，然后创建一个GeoDataFrame:

   from scipy.interpolate import griddata
   x, y = np.meshgrid(np.arange(transform[2], transform[2] + src.width * transform[0], transform[0]),
                      np.arange(transform[5], transform[5] - src.height * transform[4], -transform[4]))
   interpolated_points = griddata((x.flatten(), y.flatten()), dem_data.flatten(), (x, y), method='cubic')
   elevation_df = gpd.GeoDataFrame({'elevation': interpolated_points.reshape(x.shape),
                                  'geometry': gpd.points_from_xy(x, y)})

6. **绘图显示**:
最后，你可以使用GeoPandas绘制等值线图，比如用`matplotlib`：

   ax = elevation_df.plot(column='elevation', cmap='terrain', legend=True)
   ax.set_title("1000m等高线")