使用 RPC 进行正射校正

该存储库包含一些方法和实用程序，可帮助对原始卫星图像（例如 worldview2/3/4(RIP) 全色图像）执行正射校正。

这些算法在 2020 年 2 月 27 日的夏洛茨维尔数据科学聚会上进行了深入介绍。原始演示材料包含在/presentation

请注意，所有这些功能已经存在于 GDAL 等库中。这个代码库的目标是展示和深入研究这些子程序。

几乎所有代码都支持 Numba JIT 以获得最佳性能。还有一个正射校正例程的纯 C++ 实现，可以选择使用cppyy.

功能

ortho_tools.py

def unpack_rpc_parameters(dataset: gdal.Dataset) -> RPCCoeffs:

如果该数据集包含 RPC，则在提供 GDAL 数据集时将 RPC 系数集合作为 NamedTuple 返回

def retrieve_dem( min_lon: float, min_lat: float, degrees_lon: float, degrees_lat: float, sampling_rate: int = 1, output_path: str = "/tmp/elevation.dem", ) -> Tuple[np.ndarray, np.array]:

在给定左上角世界坐标 (min_lon, min_lat) 和以度为单位的宽度/高度的情况下加载 DEM 瓦片。

返回对象是一个元组，分别包含 DEM 图像及其 GeoTransform。

def lon_lat_to_pixel(lon: float, lat: float, geot: np.array) -> Tuple[float, float]:

使用 GeoTransform 将世界坐标重新投影到像素空间中，（6 个浮点数的数组）这对于查询 DEM 信息或查找已正射校正的图像的图像像素坐标很有用。

def linear_interp(x: float, y: float, source: np.ndarray, source_height: int) -> int: 给定像素空间坐标 x,y 和作为扁平数组的源图像（使用步幅！），返回双线性插值测量。适用于 DEM 插值以及正射校正过程中的插值。

def lon_lat_alt_to_xy( lon: float, lat: float, alt: float, rpcs: RPCCoeffs, ) -> Tuple[float, float]:

使用提供的 RPC 系数返回对应于提供的世界坐标 (lon, lat, alt) 的图像像素坐标 (x, y)

def make_ortho( x1: float, x2: float, y1: float, y2: float, width: int, source: np.ndarray, rpcs: RPCCoeffs, dem: np.ndarray, dem_geot: np.array, ) -> Tuple[np.array, float, float, float]: 使用提供的边界框 (x1, x2, y1, y2) 创建正交，即min_lon, max_lon, min_lat, max_lat, 所需宽度（像素数）、源图像、源图像的 RPC、DEM 和 DEM 的仿射 GeoTransform .

缩放工具.py

def fracture_polygon_north_up(poly: Polygon, factor_x: int, factor_y) -> Sequence[Polygon]: 将 Shapely 多边形分解为小北向上的矩形，尝试用较小的正方形填充尽可能多的区域。这些因素控制每个维度的矩形数量（分辨率）

def fracture_parallelogram(poly: Polygon, factor: int) -> Sequence[Polygon]:

将平行四边形断裂成具有非规则方向的较小平行四边形。每个单独的小平行四边形相对于原始平行四边形的方向保持不变。

def rescale_elevation_data(elevation_data: np.ndarray) -> np.ndarray: 将高程数据映射为适合查看的 0-255 8 位表示

def reproject_with_affine( coords: Sequence[Tuple[float, float]], geo_transform: np.array, resize_factor: float = 1.0, ) -> Sequence[Tuple[float, float]]: 与 基本相同ortho_tools.lon_lat_to_pixel，但适用于一系列坐标，例如，这可能正是您从 ashapely.Polygon.exterior.coords中得到的。您可以选择传入 aresize_factor以根据需要重新调整结果。

def overlay_polygon( img: np.ndarray, polygon: Polygon, color: Tuple[int, int, int, int] = (0, 255, 0, 0), opacity: float = 0.2 ) -> np.ndarray:

使用 OpenCV 将多边形刻录到图像上。例子：

def generate_triangle_mesh(elevation_data: np.ndarray, reach: int = 10) -> Tuple[np.array, np.array, np.array]: 使用非常简单且非最佳但非常快速的过程从高程数据生成三角形网格生成 3D 顶点列表作为 np.array 的元组。结果数据可以使用 OpenCV 或类似的东西进行可视化。

io_tools.py

def save_raster_as_geotiff(ortho: np.ndarray, ul_lon: float, ul_lat: float, gsd: float, filename: str) -> None: 您提供一张图像，它的左上角世界坐标（ul_lon，ul_lat）和它的 GSD（地面采样距离）（以度为单位），它会写出一个映射兼容的 GeoTIFF 文件。这可以导入 QGIS/ArcGIS 或其他任何东西以在地图上可视化正射积。

缩放工具.py

def resize(img: np.array, factor: int) -> np.array: OpenCV 调整大小

def convert_to_8bit(img: np.array) -> np.array: 将更高位的图像像素卡在 0-255 范围内……不优雅。

def gaussian_rescale(img: np.array, bitness=11, stdev_bound=3) -> np.array: 将 >8 位图像重新缩放到 0-255 范围的更优雅的方式。必须提供图像的位数。WV3 通常是 11 或 12 位。

查看之前和之后的示例：

实验.py

仅当您很勇敢时才使用它。Requirescppyy安装起来很麻烦，所以它没有被标记为必需的包。

def ortho_cpp( x1: float, x2: float, y1: float, y2: float, width: int, rpcs: RPCCoeffs, source: np.ndarray, dem: np.ndarray, dem_geot: np.array, ):

这是与 完全相同的方法签名ortho_tools.make_ortho，但它使用纯 C++ 代码作为加速版本。Numba通常大约是JITted 版本的两倍ortho_tools.make_ortho

即将推出experimental.py：使用 CUDA 正射校正您的 RPC 图像！