知识 | 常用的（python）遥感工具

介绍  

目前估计有906 颗地球观测卫星在轨运行，每天为科学和工业提供数 TB 的数据。这些卫星与雷达和光学传感器一起运行，覆盖不同的光谱范围，具有不同的光谱、空间和时间分辨率。由于广泛的地理空间数据，有可能在许多工业和政府机构中找到遥感方法的新应用。

著名的卫星系统和程序包括 ESA 的 Sentinel-1（雷达）和 Sentinel-2（光学）、NASA 的 Landsat（光学）、高分1号（光学）、高分1号（光学）、高分3号（雷达）、高分4号（光学）、高分6号（光学）、资源系列卫星（光学）等等。商业卫星不一一列举

基本上有两种类型的地理空间数据：栅格数据和矢量数据。

栅格数据  

栅格数据是由规则间隔的像素组成的网格，其中每个像素都与一个地理位置相关联，并表示为一个矩阵。像素值取决于存储的信息类型，例如数字图像的亮度值或热图像的温度值。像素的大小也决定了光栅的空间分辨率。因此，地理空间栅格数据用于表示卫星图像。光栅图像通常包含多个波段或通道，例如红色、绿色和蓝色通道。在卫星数据中，通常也有红外线和/或紫外线波段。

矢量数据  

矢量数据表示地球表面的地理特征，如城市、国界、道路、水体、产权等。这些特征由一个或多个相连的顶点表示，其中一个顶点通过x定义空间中的一个位置-、y 和 z 值。单个顶点是一个点，多个相连的顶点是一条线，多个（>3）相连且闭合的顶点称为多边形。x、y 和 z 值始终与作为元信息存储在矢量文件中的相应坐标参考系 (CRS) 相关。矢量数据最常见的文件格式是 GeoJSON、KML 和 SHAPEFILE。

为了处理和分析这些数据，需要各种工具。在下文中，我想介绍经常使用的工具。我将这些工具一一介绍，分为以下几个部分：

处理栅格数据：

•Rasterio

•Pyproj

•SNAP

•pyroSAR

•Rioxarray

           

处理矢量数据：

•Shapely

•Python-geojson

•Geojson.io

•Geopandas

•Fiona

gis开源项目：

•QGIS

•GeoServer

•Leafmap

处理气象卫星数据：

•Wetterdienst

•Wradlib

处理栅格数据  

Rasterio是一个很多模块是基于GDAL的 Python 包，可用于处理地理空间栅格数据，例如 GeoTIFF文件。为此，可以使用许多模块和函数，例如，处理来自卫星的原始数据、读取栅格数据、检索地理元数据、转换坐标、裁剪图像、合并多个图像以及以其他格式保存数据。大量的功能和易于实施使 Rasterio 成为卫星数据分析的标准工具。

建议去读下Rasterio的源代码，收货会满满的。

下面的代码示例展示了如何使用 Rasterio 读取 Sentinel-2 场景并将其显示并保存为 RGB 图像：

import numpy as np             import rasterio             from rasterio import plot                          file = 'S2BMSIL2A20210903T144719N0301R139T19LDG20210903T184347.tif'             # 打开tif文件，读取红绿蓝三个波段             with rasterio.open(file) as src:                        b, g, r = src.read(1), src.read(2), src.read(3)                        rgb = np.array([r, g, b])                        meta = src.meta.copy()                          # 更新             meta.update({'count': 3})                          # 输出新的tif             with rasterio.open(file[:-4] + '_rgb.tif', 'w', **meta) as dest:                        dest.write(rgb)                          # 可视化             rgb[rgb > 3000] = 3000             img = ((rgb / rgb.max()) * 255).astype(np.uint8)             plot.show(img)

输出应如下所示：

           

           

Pyproj  

Python 库pyproj是 PROJ 的 Python 接口，PROJ 是一个将地理坐标从一个坐标参考系统 (CRS) 转换到另一个坐标参考系统的软件库。支持大地测量变换和制图投影。Pyproj 是一个非常简单的工具，在遥感项目中必不可少，可以统一可能由所用不同数据源产生的坐标参考系统。

以下示例显示了如何对上面的 Sentinel-2 图像的边界执行简单的坐标变换：

import rasterio               import pyproj                              file = 'S2B_MSIL2A_20210903T144719_N0301_R139_T19LDG_20210903T184347.tif'               # 读取影像               # Sentinel-2地理编码代码是 EPS:32719               with rasterio.open(file) as src:                   bounds = src.bounds                   crs_from = src.crs                              # 定义要转换的坐标系 4326是WGS84坐标系               epsg_id = 4326               crs_to = pyproj.CRS.from_epsg(epsg_id)                              # 转换               proj = pyproj.Transformer.from_crs(crs_from=crs_from, crs_to=crs_to)               bounds_4326 = proj.transform_bounds(*bounds)                              print('bounds in EPSG 32719:', bounds)               print('bounds in EPSG 4326:', bounds_4326)

输出是：

bounds in EPSG 32719: BoundingBox(left=399960.0, bottom=8590240.0, right=509760.0, top=8700040.0)                bounds in EPSG 4326: (-12.752446747147742, -69.92157939546101, -11.7580241630838, -68.91008576435853)

SNAP  

SNAP全称Sentinel Application Platform，是欧空局开发的卫星数据科学探索通用工具箱。SNAP 为处理、建模和可视化卫星图像提供了一个直观的平台，特别是对于哨兵任务。该软件经过优化，可以处理大量卫星数据。它还有助于处理 SAR 数据，例如来自 Sentinel-1 的数据。

SNAP 软件安装了Sen2Cor工具。Sen2Cor 是一个为 Sentinel-2 图像执行大气校正的处理器，即将 1C 级大气顶层 (TOA) 产品转换为 2A 级大气底部 (BOA) 反射率产品。这些产品的技术规格可在此处找到：Level-1C和Level-2A。

snap是开源的，语言是纯JAVA编写的。我尝试过去编译SNAP，但是总是失败，卡在下载它的第三方包的环节，总之是莫名的bug阻止了我去看它的源代码、运行调用、逻辑结构等等。

让我心心念念的是它的S1工具箱，也就是专门去处理SAR数据和生成INSAR的工具箱。

pyroSAR   

不过不要紧，除了SNAP之外，还可以使用 Python 库pyroSAR可以处理合成孔径雷达 (SAR) 卫星数据，例如 Sentinel-1、TerraSAR-X、TanDEM-X 或 ALOS-2/PALSAR-2。访问专为 Sentinel-1 开发的 SNAP 和 GAMMA 遥感工具和功能的方法可用于此目的。

Rioxarray  

xarray是一个 Python 包，它为数组提供标签，例如尺寸、坐标和其他特定属性。因此，它使大维数组的工作更加直观。Rioxarray结合了 rasterio 的功能和 xarray 的所有优点。

Rioxarray 简化了大尺寸光栅文件的工作。它很好地概述了光栅文件的属性，例如边界、尺寸、坐标参考系统等。不幸的是，这个包没有很好的记录。而对于卫星数据的更复杂的操作，仍然需要使用栅格。

处理矢量数据  

Shapely  

Shapely是一个用于分析和操作平面地理空间特征的 Python 库。不同的地理对象由单独的类表示。点由 Point 类定义，曲线由 LineString 和 LinearRing 类定义，曲面由 Polygon 类定义。对象集合由 MultiPoint、MultiLineString 和 MultiPolygon 类实现。这些类是处理和操作矢量数据的标准，也被其他软件库和 Python 包使用。

Shapely 具有多种内置方法，可以确定对象的几何属性（例如面积、边界和长度），还可以确定与其他对象的关系，例如一个对象是否包含另一个对象、一个对象是否与另一个对象相交、是否一个物体接触另一个物体，等等。Shapely 本身没有读取矢量数据的功能，为此需要其他库，如 python-geojson、fiona 或 geopandas（见下文）。然而，将地理空间特征转换为形状物体很简单。

在下面的代码示例中，使用 GeoJSON 文件中的功能演示了上面介绍的一些标准方法：

import shapely.geometry               import geojson                              file = 'natural_areas.geojson'               # 读取 GeoJSON 文件               with open('natural_areas.geojson') as f:                   natural_areas = geojson.load(f)                              natural_areas_shapely = []               # 遍历每个图层对象               for polygon in natural_areas['features']:                   p = shapely.geometry.asShape(polygon['geometry'])                                      if p.is_valid:                     natural_areas_shapely.append(p)                   else:                     natural_areas_shapely.append(p.buffer(0))                              print('数据类型:', natural_areas_shapely[1].geom_type)               print('面积:', natural_areas_shapely[1].area)               print('边界:', natural_areas_shapely[1].bounds)               print('是否包含点', natural_areas_shapely[1].contains(shapely.geometry.Point(-76.86, -6.81)))

输出结果是：

数据类型: Polygon                面积: 0.010795427253500077                边界: (-76.952331, -6.932597, -76.794716, -6.783479)                是否包含点: True

Python-geojson  

GeoJSON 是一种标准格式，用于存储地理空间特征，如点、线或多边形及其坐标。同名的Python包提供了读取此数据、处理它并将其保存为 GeoJSON 的函数。由于其易用性，geojson 在处理以 GeoJSON 形式提供的矢量数据时是一个非常有用的工具。

上面的代码示例演示了如何使用 geojson 包读取 GeoJSON 文件。

Geopandas  

Geopandas是 pandas 的扩展，是用于数据分析和操作的标准库。pandas 已知的数据类型和方法通过基于 shapely 的地理和几何对象和方法进行了扩展。Geopandas 还具有读取和存储矢量数据的功能。这些功能取决于 fiona（见下文）。

在 geopandas 中，有形状的对象可以分组在 GeoSeries 中，可以应用各种几何方法，也可以应用经典操作进行数据分析和操作。此外，geopandas 有实用的地理编码工具。

Fiona  

使用fiona ，可以读取GDAL支持的所有光栅和矢量文件格式，也可以以此格式保存数据。除此之外，fiona 还具有一些用于几何分析的基本功能，例如确定对象的边界，以及依赖于 shapely 和 pyproj 的坐标转换。

gis开源项目  

QGIS  

QGIS是一个免费的桌面地理信息系统。它是允许显示、编辑、创建、分析和发布地理信息（栅格和矢量数据）的软件。QGIS是基于QT编写的一个GIS软件，它还提供了很多python二次开发的例子。

GeoServer  

GeoServer是一项免费服务，可用于发布地理数据。它使用开放地理空间联盟标准协议网络地图服务 (WMS)、网络要素服务 (WFS)、网络覆盖服务 (WCS)、网络处理服务 (WPC) 和网络地图瓦片服务 (WMTS)。

             

地理空间栅格和矢量数据可以使用上述协议发布，随后用于 Web GIS 应用程序（例如通过 Leaflet、Mapbox 或 OpenLayers）或 GIS 软件（例如 QGIS 或 ArcGIS）。这能够以交互方式提供地理数据，而无需长时间加载和大量数据传输。

Leafmap  

Leafmap是一个 Python 包，用于在 Jupyter 笔记本中对地理空间数据进行交互式可视化。只需几行代码，就可以轻松设计交互式地图。Leafmap 非常强大，因为它构建并集成了各种地理空间映射和分析工具，例如ipyleaflet、whiteboxtools和folium。

Leafmap 支持所有常见的矢量和栅格数据格式。您可以从不同的基础地图中进行选择，包括卫星地图、混合地图、街道、地形等等。可以绘制自己的多边形并将它们导出为矢量文件或将它们用于后续操作。该软件包还提供上传功能。这些地图是高度可配置的，因为多边形、线条、标记、弹出窗口和图例存在许多样式选项。如果您将地理空间数据托管在图块服务器上，例如通过 Geoserver，leafmap 可以从 WMS 或 WFS 等服务流式传输栅格和矢量数据。

             

吴秋生老师写的GEEMAP使用了Leafmap包。

             

GEE  

GEE，google earth engine。是遥感界的大杀器。提供了JS和PYTHON两个API接口，实现了遥感大数据可视化，对于研究一些地球物质长时间序列变化的课题提供了巨大帮助。

处理气象卫星数据  

Wetterdienst  

Wetterdienst是一个 Python 包，可用于从德国气象局 (DWD) 检索气象数据。其中包括历史天气数据、基于天气模型的预报和雷达数据（RADOLAN 和 RADVOR 等）。简单雷达数据（如温度、大气降水等）的时间分辨率从月到年不等。RADOLAN 数据以分钟到每日的时间分辨率提供。

不幸的是，目前没有内置函数来创建 RADOLAN 数据的时间序列，而这些数据通常用于分析天气数据。

Wradlib  

Wradlib是一个 Python 库，可用于处理各种天气雷达数据，包括 RADOLAN、DX 和 BUFR。特别有用的是天气数据可以转换为合适的栅格格式，例如 GeoTIFF。

             

完结  

由于篇幅原因，本文没有将所有遥感相关的PYTHON工具介绍。

关注我，分享更多有趣、好玩、有用的东西。

个人公众号：remote sensing