引言

2020年8月中旬左右，由于700多次云对地雷击，加州中部和北部发生了一系列野火(350+)。美国国家航空航天局™的索米核电站卫星在这些地点上空拍摄到了引人注目的森林大火夜间画面。卫星在探测野火方面非常有用，不仅如此，它们还跟踪烟雾的传播，并有助于确定破坏的程度。NASA’s Suomi NPP Satellite

每天都有数十亿项活动，世界日益相互联系，我们正以前所未有的速度发展。卫星图像将在跟踪这些全球变化方面发挥关键作用。卫星数据的应用仅限于我们的想象。这些数据正在被私人和政府机构大规模使用。监测灾害，跟踪农业发展，分析港口的航运活动，观察城市化的增长，这些都是很少的应用。

与标准的RGB图像相比，卫星图像有很多优点。它不仅限于可见光光谱，还可以使用其他电磁波，如红外线，来创建RGB以外的波段。此外，即使在夜间或恶劣天气下，也可以进行成像。这些图像通常以Geotiff文件格式存储。此格式允许将地理空间信息与图像数据一起嵌入到文件中。在本博客中，我们将通过处理为此问题陈述链接提供的数据集中的图像来了解如何使用此文件格式。目标是检测提供的卫星数据中的建筑物足迹。在探索其中一幅卫星图像的同时，我们将致力于回答以下问题：link

• 如何浏览嵌入在tif文件中的不同地理空间信息？
• 如何将建筑物边界从一个CRS转换到另一个CRS？
• 如何使用Geojson文件中提供的建筑物外形从给定的图像生成掩模？
• 如何准备以非常规格式TIF存储的卫星图像，以供消费到AI模型(欧元)中，从图像中创建瓷砖？
• 如何获得原始图像的较低分辨率版本？

概念概述

输入：

• 图像的TIF文件和包含建筑物迹线坐标的Geojson文件。
• TIF包含图像数据以及存储在配置文件属性中的图像的元数据
• 本博客中使用的几个示例文件可以在此处的链接中找到

杰森：

• 具有三个主键的JSON文件-EUROPE、CRS和FEATURES
• 要素键包含要素列表，其中每个要素都是一个字典，其中包含建筑物轮廓线的信息(如轮廓线坐标)的™
• CRS键提到存储轮廓线坐标时使用的CRS

使用的库：

由于TIF图像是大文件，并且包含重要的元数据，因此我们将使用专门为处理和操作这些地理空间数据而设计的库。

• Rasterio：读/写TIF，访问配置文件和其他元数据。创建切片并将其存储为TIF
• AERONET：高效地读取、解析和存储Geojson中的信息。将批注从一个CRS重新投影到另一个CRS

• 地理坐标系：横跨整个地球的坐标系(例如纬度/经度)。
• 投影坐标系：为使特定区域的视觉失真最小化而定位的坐标系(例如，Robinson、UTM、州平面)

为什么选择多个CR？

CRS经过优化，最能代表以下各项：

• 形状和/或
• 缩放距离和/或
• 面积

CRS要么优化形状，要么优化距离，要么优化面积，没有一个CRS在优化所有这三个元素方面做得很好。一些CRSâeuro™还针对特定地区(欧元)进行了优化，例如美国或欧洲。

有许多格式可用于记录CRS。Proj4、WKT(熟知文本)和EPSG。EPSG得到了广泛的应用。下面我们来了解一下。

EPSG简介

• EPSG代码是指以EPSG提供的格式描述的不同CR
• 每个CRS都定义了它使用的原点和它选择用来描述地球的地球形状，还定义了它将用来引用该参考系上的坐标的单位制
• 例如，WDS84或WGS 1984使用赤道和本初子午线的交点作为原点，并使用度(经度)作为单位来描述坐标
• EPSG：4326是一种未投影的地理CRS，而下面使用的ESPG：32327是一种投影的CRS。
• 顾名思义，未投影的CRS将地球视为3D，并不将其投影到2D表面上，而投影的CRS则使用3D到2D的映射系统。可以通过多种方式完成3D到2D的转换。并且没有办法完美地映射它，因此在区域、形状和比例链接上总是会有一些失真(或折衷
• 这里是将一条坐标从一条CRS转换到另一条链路链接

Rasterioâuro“处理卫星图像。

import os
import json
import random
import rasterio
import numpy as np
import aeronet.dataset as ds
import matplotlib.pyplot as plt


from PIL import Image
from glob import glob
from tqdm import tqdm
from fastcore.all import *
from typing import Union
from rasterio.features import geometry_mask
from rasterio.windows import Window
from rasterio.warp import calculate_default_transform, reproject, Resampling, transform_geom

我们使用的是来自Spacenet数据集的示例图像

src_path = "buildseg/train_tier_1/train_tier_1/znz/33cae6/33cae6.tif"

浏览TIF文件中存在的地理空间信息

> data = rasterio.open(src_path)
> type(data)
[Out] rasterio.io.DatasetReader

• Rasterio为我们提供了一个非常有用的数据加载器。我们可以很容易地访问数据中存在的各种地理空间信息。让我们来看看重要的几个

配置文件包含如下信息：

• 仿射变换：帮助将CRS转换为像素，并将像素转换为CRS
• CRS：以像素为单位的坐标参考系宽度和高度
• 块大小和块大小和平铺：用于定义内部平铺。一些我们将要实施的事情。
• nodata：数据中是否缺少像素
• 计数：波段数

在创建TIF图像时，提供配置文件和图像数据很重要。这一点很重要，在创建切片时将重新访问它

> data.profile
[Out] {'driver': 'GTiff', 'dtype': 'uint8', 'nodata': None, 'width': 37113, 'height': 34306, 'count': 4, 'crs': CRS.from_epsg(32737), 'transform': Affine(0.0774800032377243, 0.0, 531847.0,
       0.0, -0.0774800032377243, 9367848.0), 'blockxsize': 512, 'blockysize': 512, 'tiled': True, 'compress': 'jpeg', 'interleave': 'pixel'}

• 边界给出了给定CRS系统中图像的左下角、最左上角和最右上角的坐标。此图使用epsg 32737编码，单位为米。

> data.bounds
[Out] BoundingBox(left=531847.0, bottom=9365189.971008927, right=534722.5153601617, top=9367848.0)

• 我们还可以得到图像中心的长纬度。

> data.lnglat()
[Out] (39.30061294996244, -5.731031019254337)

• 我们还可以获得以米/像素为单位的垂直和水平分辨率

> data.res
[Out] (0.0774800032377243, 0.0774800032377243)

降低图像的分辨率

• 原始图像的分辨率很高，降低分辨率会加快处理速度。进行分析不需要极高的分辨率，因为它会减慢处理步骤。

现在我们将创建分辨率为0.1的等效tif

> dst_res = 0.1

• 在下面的所有参数中-欧元“压缩”、平铺是很重要的

>DST_PATH=“随机2.tif”

> command = f"gdalwarp -co COMPRESS=JPEG -co TILED=YES  -co NUM_THREADS=ALL_CPUS -r bilinear " + \
                   f"-tr {dst_res} -{dst_res} {src_path} {dst_path}"
> os.system(command)  

• 分辨率降低了，所以现在每个像素覆盖了更多的区域，因此我们注意到高度和宽度都减小了

> img = rasterio.open(dst_path)
> img.profile
[Out] {'driver': 'GTiff', 'dtype': 'uint8', 'nodata': None, 'width': 28755, 'height': 26580, 'count': 4, 'crs': CRS.from_epsg(32737), 'transform': Affine(0.1, 0.0, 531847.0,
       0.0, -0.1, 9367848.0), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'compress': 'jpeg', 'interleave': 'pixel'}

可视化TIF图像

• 我们被提供了四个乐队。我们会选择RGB波段并将其可视化。这些图像看起来很壮观。我们可以很容易地将图像切片并仔细观察，但稍后我们将通过将图像分割为瓷砖来实现这一点。

plt.figure(figsize=(22, 22))
plt.imshow(img.read([1,2,3]).transpose(1, 2, 0))

探索Geojson

• 我们将探索相应tif文件的Geojson
• 我们将修复图像和Geojson之间的epsg代码中的不匹配问题。
• 使用Geojson中存在的数据，我们将使用降低的分辨率创建新图像的掩膜

base_dir = os.path.dirname(src_path)
di, name = os.path.split(base_dir)
geojson_dir_name = name + "-labels"
geojson_file_name = name + ".geojson"

FeatureCollection是Aeronet.DataSet模块的一种方法

• 专门为阅读和处理Geojson而制作
• FC是Geojson中存在的所有要素(建筑物轮廓线)的矢量
• 它还提供了我们将在下面讨论的其他属性
• 可以处理CRS转换的Geojson数据加载器
• 还记得我们讨论过有多个CR吗？现在，我们需要确保我们的图像坐标和Geojson中的坐标在相同的CRS中表示。如果没有，我们将无法映射它们并创建蒙版
• 遗憾的是，默认的dataloader FeatureCollection.read中硬编码了CRS，因此在读取Geojson时不准确
• 此外，为了简化Geojson从一个CRS到另一个CRS的转换，我们对原始数据加载器进行了一些更改，并修补了一个新的FeatureCollection.read_Generic
• patch_to允许我们将函数修补到给定的类，如下所示

代码的解释超出了本博客的范围，我们将在下一篇博客中讨论这一点

> fc = ds.FeatureCollection.read_generic(geo_path)
[Out] urn:ogc:def:crs:OGC:1.3:CRS84

让我们对照图像中使用的CRS引用来查看Geojson的默认CRS

> fc.crs, img.crs
[Out] ('urn:ogc:def:crs:OGC:1.3:CRS84', CRS.from_epsg(32737))

修复epsg不匹配并将坐标从一个CRS转换为另一个CRS

• 拥有不同的CRS的问题是，我们将有不同的框架来表示地球上的同一个点
• 就像EPSG：32737一样，基本单位是米，而EPSG：4326度。而且两者都有不同的起源参照
• 由于米很容易解释，我们将Geojson中的坐标转换为epsg：32737。
• fc.reproject提供了一种将其Geojson从一个CRS映射到另一个CRS的方法。这是上天赐予的礼物™。
• 原始Geojson不会被覆盖，但会处理存储在矢量中的每个要素，并且在打印每个要素的几何时可以观察到更改，如下所示

> fc = fc.reproject(img.crs)
> fc.crs
[Out] CRS.from_epsg(32737)

更简单的方法是使用新的dataloader加载Geojson，然后执行以下操作

> fc = ds.FeatureCollection.read_generic(geo_path, dst_crs=img.crs)
[Out] urn:ogc:def:crs:OGC:1.3:CRS84

你会得到这个

> fc.crs
[Out] CRS.from_epsg(32737)

• 每个要素表示单个建筑物覆盖区
• fc.Feature从所有掩码(建筑物覆盖区)中提取？euroœproperties？euro�密钥

> fc.features[0]

让我们探索第一个功能(建筑物占地面积)。

图像中的建筑物覆盖区总数为

> len(fc)
4439

• 几何包含建筑物迹线的面的所有坐标

> fc[0].geometry
[Out] {'type': 'Polygon',
 'coordinates': (((533700.1052929291, 9367442.727510184),
   (533702.2944919249, 9367446.207351241),
   (533704.6579756252, 9367445.290231245),
   (533707.6459179376, 9367450.135843525),
   (533710.2917512723, 9367450.400426859),
   (533713.5990429427, 9367449.143656025),
   (533715.2526887775, 9367446.828551855),
   (533755.0559440283, 9367429.762926837),
   (533754.857506527, 9367428.836885171),
   (533770.1041211254, 9367420.00641641),
   (533766.6645377881, 9367418.352770574),
   (533761.3728711164, 9367421.593916412),
   (533753.369225275, 9367421.263187245),
   (533713.4667512759, 9367440.908499764),
   (533708.654641896, 9367433.814359132),
   (533702.30464189, 9367437.386234136),
   (533703.3629752239, 9367441.222692475),
   (533700.1052929291, 9367442.727510184)),)}

• 此建筑迹线覆盖的面积

> fc[0].shape.area
[Out] 531.2663747649708

• 提供的元数据

fc[0].properties

从Geojson中显示的几何创建遮罩

GEOMETRY_MASK：使用Geojson文件中存在的几何图形列表创建整个图像的蒙版。我们需要提供我们正在为其创建蒙版的图像的变换，因为它内部组合了所有几何图形并将其映射出来

> geometries = [f.geometry for f in fc]

• 通过设置掩码的形状，并通过对图像的变换，得到数字数组形式的掩码

> profile = new_image_src.profile.copy()
> profile

将蒙版存储为tif文件并将其可视化

> mask_path = "mask.tif"
> type(mask)
[Out] numpy.ndarray

• 我们将掩码数据存储为tif，这样我们还可以对地理空间信息进行编码。因此，我们需要为其创建配置文件。
• 纵断面与光栅数据中显示的方法相同，如下所示。
• 一个重要的关键是“EUROUœTransform”EURO�。您的配置文件需要有转换密钥。如果没有它，您就不能将图像保存为栅格数据，™t/不应该将™另存为栅格数据。

> img.profile
[Out] {'driver': 'GTiff', 'dtype': 'uint8', 'nodata': None, 'width': 28755, 'height': 26580, 'count': 4, 'crs': CRS.from_epsg(32737), 'transform': Affine(0.1, 0.0, 531847.0,
       0.0, -0.1, 9367848.0), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'compress': 'jpeg', 'interleave': 'pixel'}

• 创建蒙版的配置文件
• 稍后我们将在创建切片时重新访问变换部分

> profile = dict(
        driver="GTiff",
        height=mask.shape[0],
        width=mask.shape[1],
        count=1,
        crs=profile["crs"],
        transform=profile["transform"],
        dtype=mask.dtype,
        NBITS=1,
    )

让我们把面具形象化

> plt.figure(figsize=(22, 22))
> plt.imshow(rasterio.open(mask_path).read()[0])

为图像和蒙版创建平铺

• 对于所有的深度学习专家来说，在你头脑正常的情况下，你永远不会将30kx30k大小的图像输入到模型中
• 因此，我们将把图像拆分成小块，并仔细地为小块生成蒙版
• 我们将切片存储为tif，因此需要为每个切片设置配置文件
• 为了提取每个平铺的像素值，我们将使用Rasterio中的窗口属性。不清楚吗？呆在原地别动。
• 用于存储平铺图像和蒙版的目录

dst_image_dir = "images"
dst_mask_dir =  "masks"

平铺图像

将像素坐标转换为CRS坐标

为什么我们需要了解这种转换以生成平铺？

• 还记得您需要配置文件dict将Numpy数组存储为栅格吗？是的，这就是我们在这里需要它的原因。
• 我们可以使用滑动窗口策略轻松地对原始图像的模糊数组进行切片，但我们也希望将切片存储为栅格
• 因此，我们需要创建所有瓷砖的配置文件字典。这意味着我们还需要创建变换矩阵。
• 让我们通过理解转换来理解这个转换矩阵
• 我们将使用父图像的变换矩阵为每个切片创建变换矩阵
• 在为每个平铺创建变换矩阵时，我们需要使用提供的分辨率计算给定CRS系统中的最左上角坐标
• 对于每个切片，如果将原点固定到最左上点，则可以计算沿两个轴到原点的距离为
• X轴：X*转换.a
• Y轴：Y*trans.e
• x和y分别是沿x和y轴遍历的像素数。
• 从左向右遍历被视为正向平移。从北向南穿越被认为是负翻译。
• 一旦我们以米为单位遍历了距离，我们将把原点移回父图像的最左上角，以计算坐标
• 对于计算最左边的coord，我们将使用：coord_x=trans.c+x*trans.a
• 为了计算最上面的coord，我们将使用公式eq：coord_y=trans.f+y*trans.e
• 最上面的坐标形成平铺变换矩阵的c&e元素

img.profile

平铺的大小=1024

size = 1024

• 以下内容将有助于命名切片

i = 0
city = src_path.split("/")[-3]
code = src_path.split("/")[-2]
print(city, code)

使用滑动窗口(带跨距)策略并存储切片

• 这里SIZE=STRIDE，但是如果我们想要改变大小，那么我们只需在两个for循环中用range的step参数中所需的步长替换SIZE
• 栅格的.read属性中的无边界参数有助于处理窗口将部分延伸到图像之外的角点情况

Similarly we will split the mask tif

可视化随机平铺和蒙版

for index in random.sample(range(754), 5):
    name_index = "{}_{}_{}.tif".format(city, name,     str(index).zfill(5))
    img_tile = rasterio.open(os.path.join(dst_image_dir, name_index))
    mask_tile = rasterio.open(os.path.join(dst_mask_dir, name_index))
    plt.figure(figsize=(12,12))
    plt.subplot(1,2,1)
    plt.imshow(img_tile.read().transpose(1, 2, 0))
    plt.subplot(1,2,2)
    plt.imshow(mask_tile.read()[0])

声望与深度学习

• 一旦我们有了平铺图像和蒙版，我们就可以训练Mask-rcnn并构建构建分割模型
• 此检测器将以瓷砖为输入，并检测建筑物占地面积
• 我们可以对测试tif文件的平铺运行推断，并重新缝合它，以通过检测来可视化整个mosiac。

作者：？EURO“库纳尔·辛格和普拉卡什·杰伊

资源