Analyzing Geography Data (Beginners' Tutorial)

数据来源

卫星类型

  1. Sentinel-2:提供混合分辨率的13 Bands MSI。分辨率有$10m \times 10m$,$20m \times 20m$,$60m \times 60m$,bands中心波长从442.3nm到2185.7nm。

    image-20201213154533413
  2. Planet:这个数据源提供$3m \times 3m$的4 bands MSI image。分别是R、G、B、NIR。

    • 这个数据源似乎只对美国国内的院校机构开放,不过不是很确定,需要的可以试试。

    Planet Specification

  3. RapidEyes

Sentinel-2 Data Source

  1. USGS Earth Explorer: Link. It support downloading data within five years. And it not only support sentinel-2 data, but also contains many other satellite datasource.

    image-20201227143234319

  2. Copernicus Open Access Hub: Link. It supports only one year’s old data. You can request the older data, but not guarante the fetch time.

    image-20201227145324164

  3. Amazon AWS Sentinel-2 Service: Link

Planet Data Source

Specification

  1. Sentinel-2 MultiSpectral Instrument (MSI) Overview
  2. Sentinel-2 Specification Doc
  3. Planet Specification

软件与Python包简介

  • GDAL –> Fundamental package for processing vector and raster data formats (many modules below depend on this). Used for raster processing.
  • Rasterio –> Clean and fast and geospatial raster I/O for Python. Guidebook.
  • Geopandas –> Working with geospatial data in Python made easier, combines the capabilities of pandas and shapely.
  • Shapely –> Python package for manipulation and analysis of planar geometric objects (based on widely deployed GEOS).
  • Fiona –> Reading and writing spatial data (alternative for geopandas).
  • Pyproj –> Performs cartographic transformations and geodetic computations (based on PROJ.4).
  • Pysal –> Library of spatial analysis functions written in Python.
  • Geopy –> Geocoding library: coordinates to address <-> address to coordinates.
  • GeoViews –> Interactive Maps for the web.
  • Networkx –> Network analysis and routing in Python (e.g. Dijkstra and A* -algorithms), see this post.
  • Cartopy –> Make drawing maps for data analysis and visualisation as easy as possible.
  • Scipy.spatial –> Spatial algorithms and data structures.
  • Rtree –> Spatial indexing for Python for quick spatial lookups.
  • RSGISLib –> Remote Sensing and GIS Software Library for Python.
  • python-geojson-> Deal with geojson format files.

Sentinel-2 大气校正

Sentinel-2一般有两种standard,一个是A级别,一个是C级别。C级别的数据是你可以从任意网站上下载到的数据,它没有经过大气校正,是粗数据,每个区块的反射率可能有不同,不适合直接用作深度学习数据。

经过Sen2Cor软件校正后可以得到A级别的数据。

Sen2Cor是欧空局发布的一个软件,它即可以作为SNAP的插件安装,也可以作为独立命令行软件使用。推荐后者,更为快速、稳定,尤其适用于大量数据时,可以用脚本批处理。官网链接

查询坐标映射

  • 地球是圆的。
  • 卫星上拍的一张矩形照片所对应的区域并非是矩形的。
  • 使用卫星图片时要先将其映射到二维展开的坐标系中。
  • 整个地球被分为了许多预先订好的区域。
  • 每个区域有一个编号。
  • 编号可以在EPSG网站查询。

SNAP

  • 欧空局自己用作处理Sentinel-2的软件。Link

  • 只用来处理Sentinel-2,尤其是预处理,包括校正,reprojection,粗crop,统一各个bands分辨率等等,非常好用。因为是亲儿子,所以甚至可以直接读取Sentinel-2每个文件的压缩包,总之十分便利。

  • 比较古老,interface有年代感,功能也相较于其他软件比较局限。

  • 推荐用作第一步预处理。

  • 使用流程:读取->剪裁->correction->resize->reprojection->导出

  • 大部分需要用到的功能都在Raster->Geometric

  • 经过尝试、搜索、确认,SNAP并不提供便利的选定区域截图,或是依据shapefile剪裁,所以目前最佳的方法是,先通过zoom地图和改变窗口大小确保需要的部分大致在view的可视范围内,然后右键,选择Spatial Subset from View,然后可视区域即可被剪裁。注意,这只是粗剪裁,所以尽量多包括一点,也不要少任何一部分。剪裁后的内容可以后续在ArcGIS Pro中进一步处理。

    image-20201227155911702
  • 展示图:

    image-20201227155421761

    image-20201227155558109

ArcGIS Pro

  • 我能找到的最强大的可用于卫星图像处理的软件。

  • 专业,美观,支持format多,有着强大的raster functions。

  • 版权软件,下载之前先确认自己学校或公司是否提供License。

  • 有时导出raster会出现随机bug,导致:导出可能是纯黑的图片,导出部分没有按照期望剪裁等。很恶心,但似乎也没有更好的选择。

  • 解决导出问题:

    1. 关闭导出窗口,重新操作,多试几次。可以解决绝大部分问题。
    2. 和剪裁、mask等有关的问题可以先使用raster function进行这些操作,再直接导出前面操作的结果layer。
  • 关于剪裁后的卫星图片和原始图片有着明显亮度对比度区别的问题:

    1. 首先,这不是一个bug,而是一个feature。。。实际上的卫星图都很暗沉的,所以软件原生提供一个“显示方法”的函数,调整图片曲线,使得显示的图片比较亮,容易看清楚细节。然而,剪裁后的图片有着不一样的统计值,所以在有些显示函数下,显示的结果和剪裁前结果不同。

      image-20201227155342029
    2. 但是不用担心,因为导出时候并不会考虑这个显示函数。即:剪裁前后的导出图像数值是相同的。

  • 下面是两张展示图:

image-20201227155051608

image-20201227155259214

QGIS

  • 开源软件,支持很多插件,比如直接下载Sentinel-2数据,Sen2Cor等。Link

  • 支持很多format的数据。

  • 问题是,不稳定,效率低,容易崩溃(软件&心态),甚至左侧Explorer遇到大文件夹都要经常转很久才能进去,或是干脆就直接转崩了。我怀疑他们要分析每个文件夹所有文件之后再显示列表。。。总之,慎重。

    image-20201227160344273

Pythonic Method

虽然前面介绍了几个软件,但是说实话,处理一两个可以,批量处理几十几百甚至几十万就有点力不从心了。所以最后还是狠下心研究了一遍Python处理这些数据的方法,写出了一批适合我项目用的函数。不一定适合所有人,但可以作为参考:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
import numpy as np
import json
import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon
import rasterio as rio
from rasterio.mask import mask
from rasterio.warp import calculate_default_transform, reproject, Resampling

def bbox(shp):
""" Compute the bounding box of a certain shapefile.
"""
piece = np.array([i.bounds for i in shp['geometry']])
minx = piece[:,0].min()
miny = piece[:,1].min()
maxx = piece[:,2].max()
maxy = piece[:,3].max()
return minx, miny, maxx, maxy

def edge_length(shp):
""" Compute the x and y edge length for a ceratin shapefile.
"""
minx, miny, maxx, maxy = bbox(shp)
return round(maxx-minx,3), round(maxy-miny,3)

def shape2latlong(shp):
""" Turn the shapefile unit from meters/other units to lat/long.
"""
return shp.to_crs(epsg=4326)

def bbox_latlong(shp):
""" Compute the latitude-longitude bounding box of a certain shapefile.
"""
shp = shape2latlong(shp)
return bbox(shp)

def bbox_polygon(shp):
""" Return the rectangular Polygon bounding box of a certain shapefile.
"""
minx, miny, maxx, maxy = bbox(shp)
return Polygon([(minx, miny), (minx, maxy), (maxx,maxy), (maxx, miny)])

def merge_polygon(shp):
""" Merge a shapefile to one single polygon.
"""
return shp.dissolve(by='Id').iloc[0].geometry

def polygon2geojson(polygon):
""" Turn a polygon to a geojson format string.
This is used for rasterio mask operation.
"""
if type(polygon) == Polygon:
polygon = gpd.GeoSeries(polygon)
return [json.loads(polygon.to_json())['features'][0]['geometry']]

def sen2rgb(img, scale=30):
""" Turn the 12 channel float32 format sentinel-2 images to a RGB uint8 image.
"""
return (img[(3,2,1),]/256*scale).astype(np.uint8)

def cropbyshp(raster, shp):
""" Crop a raster using a shapefile.
"""
# Reproject the shapefile to the same crs of raster.
shp = shp.to_crs({"init": str(raster.crs)})
# Compute the rectangular Polygon bounding box of a certain shapefile.
bbpoly = bbox_polygon(shp)
# Execute the mask operation.
out_img, out_transform = mask(dataset=raster, shapes=polygon2geojson(bbpoly), crop=True, all_touched=True)
return out_img

def write_raster(raster, path):
""" Write a created raster object to file.
"""
with rio.open(
path,
'w',
**raster.meta
) as dst:
dst.write(raster.read())

def sen_reproject(src, dst_crs, out_path):
""" Reproject a raster to a new CRS coordinate, and save it in out_path.
Args:
src: Input raster.
dst_crs: Target CRS. String.
out_path: The path of the output file.
"""
transform, width, height = calculate_default_transform(
src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height
})

with rio.open(out_path, 'w', **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rio.band(src, i),
destination=rio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.cubic)

def mask_A_by_B(A, B):
""" Generate a mask from B, and applied it to A.
All 0 values are excluded.
"""
mask = B.sum(axis=0)>1e-3
masked_A = mask*A
return masked_A

Reference

  1. Python GIS 超完整教程
  2. Professor Zia’s Personal Website
  3. sentinel-2数据下载 大气校正 转ENVI格式
  4. 03-SNAP处理Sentinel-2 L2A级数据(一)
  5. UBC Course Notebook
  6. 利用Sen2cor对哨兵2号(Sentinel-2)L1C多光谱数据进行辐射定标和大气校正
  7. Documentation of Sentinel Hub Python package
  8. sentinelhub-py