Python下的GeoTrellis —— GeoPySpark

Python下的GeoTrellis —— GeoPySpark

什么是GeoPySpark

GeoPySpark是Scala库GeoTrellis的Python语言库。与GeoTrellis一样,这个项目也是在Apache 2许可下发布的。

GeoPySpark试图利用GeoTrellis来实现对栅格数据的读取、写入和操作。因此,它能够扩展到数据,并且仍然能够很好地执行。

除了栅格处理,GeoPySpark还允许将栅格数据渲染成PNG。该项目的目标之一是能够以web速度处理栅格数据并对大型数据集执行批处理。

为什么要用GeoPySpark

Python中的处理栅格已经取得了长足的进步;然而,随着数据集大小的增加,仍然会出现问题。无论是性能还是易用性,随着越来越多的数据向公众开放,这类问题将变得更加普遍。

人们可以求助于地理网格来解决上述问题(并且应该尝试一下!),但这也带来了新的挑战。Scala虽然是一门强大的语言,但它的学习曲线有些陡峭。这可能会耽误那些没有时间和/或兴趣学习一门新语言的人。

通过Scala的速度和可伸缩性以及Python的易用性,GeoPySpark可以解决这种困境。

核心概念

因为GeoPySpark是一个现有项目GeoTrellis的库,所以一些术语和数据表示形式得以延续。本节除了描述GeoPySpark中如何表示栅格类型外,还将解释这个术语。

在开始之前,本指南中的所有示例都需要导入一下包:

1
2
3
4
import datetime
import numpy as np
import geopyspark as gps
from shapely.geometry import box

栅格数据

GeoPySpark不同于其他地理空间Python库(如rasterIO)的栅格表示方式。在GeoPySpark中,它们由Tile类表示。该类包含一个numpy数组(称为cells),除了与数据相关的其他信息外,该数组还表示格网的单元格。除了cells之外,Tile还可以具有栅格数据的no_data_value

注意: GeoPySpark中的所有栅格数据都可以表示为多个波段,即使原始光栅仅包含一个波段。

1
2
3
4
5
6
7
8
9
arr = np.array([[[0, 0, 0, 0],
[1, 1, 1, 1],
[2, 2, 2, 2]]], dtype=np.int16)

# The resulting Tile will set -10 as the no_data_value for the raster
gps.Tile.from_numpy_array(numpy_array=arr, no_data_value=-10)

# The resulting Tile will have no no_data_value
gps.Tile.from_numpy_array(numpy_array=arr)

范围

描述栅格数据所代表的地球上的位置。这个区域是由坐标参考系中的坐标表示的。因此,根据所使用的坐标系统不同,范围的值可能会有所不同。Extent也可以作为一个bounding box 来引用。

注意: Extent内的值必须是float,而不是double

1
2
extent = gps.Extent(0.0, 0.0, 10.0, 10.0)
extent

坐标范围

ProjectedExtent描述了除了CRS外栅格数据在地球上所代表的区域。可以使用EPSG代码proj4字符串来指示ProjectedExtent的CRS。

1
2
3
4
5
6
# Using an EPSG code
gps.ProjectedExtent(extent=extent, epsg=3857)

# Using a Proj4 String
proj4 = "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "
gps.ProjectedExtent(extent=extent, proj4=proj4)

临时坐标范围

ProjectedExtent类似,TemporalProjectedExtent描述栅格数据所表示的地球上的区域、CRS和数据所表示的时间。这个时间点称为instant,是datetime.datetime 的一个实例。

1
2
time = datetime.datetime.now()
gps.TemporalProjectedExtent(extent=extent, instant=time, epsg=3857)

切片布局

TileLayout描述表示栅格如何在一个层中的组织和分类。layoutColslayoutRows分别详细说明了网格本身有多少列和行。而tileColstileRows 表示每个栅格具有多少列和行。

1
2
3
# Describes a layer where there are four rasters in a 2x2 grid. Each raster has 256 cols and rows.
tile_layout = gps.TileLayout(layoutCols=2, layoutRows=2, tileCols=256, tileRows=256)
tile_layout

Layout定义

LayoutDefinition描述了栅格数据在一个层中的组织方式以及栅格所覆盖的区域。

1
2
layout_definition = gps.LayoutDefinition(extent=extent, tileLayout=tile_layout)
layout_definition

切片策略

通常情况下,层的布局是未知的。这里有两种不同的瓦片生成策略,它们将根据给定的数据生成布局,而不是费力地找出最优布局。

局部策略

LocalLayout是第一种瓦片生成策略,它生成的布局是在给定瓦片大小的层内的所有像素上构建网格。生成的布局将匹配栅格内单元格的原始分辨率。

注意: 此布局不能用于创建显示层。相反,它最好用于将执行操作和分析的层。

1
2
3
4
5
6
7
8
# Creates a LocalLayout where each tile within the grid will be 256x256 pixels.
gps.LocalLayout()

# Creates a LocalLayout where each tile within the grid will be 512x512 pixels.
gps.LocalLayout(tile_size=512)

# Creates a LocalLayout where each tile within the grid will be 256x512 pixels.
gps.LocalLayout(tile_cols=256, tile_rows=512)

全局策略

另一种切片生成策略是GlobalLayout,即在全局范围CRS上构建网格。结果层的单元格分辨率乘以CRS的2次方。因此,使用这种策略将导致原始栅格数据的向上或向下采样。

注意: 这种布局策略应该使用当结果层要在TMS服务器中展现。

1
2
3
4
5
# Creates a GobalLayout instance with the default values
gps.GlobalLayout()

# Creates a GlobalLayout instance for a zoom of 12
gps.GlobalLayout(zoom=12)

可以从上面的两个示例中注意到,GlobalLayout默认情况创建布局时下不会给定的缩放级别。相反,它根据栅格内单元格的大小决定缩放的大小。如果希望为特定的缩放级别创建布局,则必须设置zoom参数。

SpatialKey

SpatialKey描述栅格在布局网格中的位置。这个网格是一个二维平面,其中栅格的位置由一对坐标colrow分别表示。顾名思义,SpatialKey只处理空间数据。

1
gps.SpatialKey(col=0, row=0)

SpaceTimeKey

SpatialKeys一样,SpaceTimeKeys描述栅格在布局中的位置。然而,网格是一个三维平面,其中光栅的位置由一对坐标colrow表示,以及一个表示时间点instant的z值表示。与TemporalProjectedExtent中的instant类似,这也是datetime.datetime的一个实例。因此,SpaceTimeKeys 处理时空数据。

1
2
time = datetime.datetime.now()
gps.SpaceTimeKey(col=0, row=0, instant=time)

界限

Bounds表示布局网格中键的范围。它同时具有minKeymaxKey属性。根据层中的数据类型,这些键可以是SpatialKeySpaceTimeKeyminKey是网格中最左边的单元格,maxKey 是网格中最右边的单元格。

1
2
3
4
5
6
7
8
9
10
11
# Creating a Bounds from SpatialKeys

min_spatial_key = gps.SpatialKey(0, 0)
max_spatial_key = gps.SpatialKey(10, 10)
bounds = gps.Bounds(min_spatial_key, max_spatial_key)
bounds

# Creating a Bounds from SpaceTimeKeys
min_space_time_key = gps.SpaceTimeKey(0, 0, 1.0)
max_space_time_key = gps.SpaceTimeKey(10, 10, 1.0)
gps.Bounds(min_space_time_key, max_space_time_key)

元数据

元数据包含层中值的信息。这些数据属于层中包含的数据的布局、投影和范围。

下面的示例展示了如何手工构造元数据,但是,这几乎从来不需要,而且可以使用更简单的方法生成元数据。对于RasterLayer,可以调用方法collect_metadata(),而TiledRasterLayer具有属性layer_metadata

1
2
3
4
5
6
7
# Creates Metadata for a layer with rasters that have a cell type of int16 with the previously defined
# bounds, crs, extent, and layout definition.
gps.Metadata(bounds=bounds,
crs=proj4,
cell_type=gps.CellType.INT16.value,
extent=extent,
layout_definition=layout_definition)

特性

特性 是带有某种关联元数据的形状几何图形。这种类型的主要目的是提供一种栅格化具有不同单元格值和可能重叠的多个几何图形的方法。

单元值

在实践中,一个Feature可以具有任何类型的元数据,但是CellValue需要与要栅格化的特性一起使用。CellValue有两个参数:valuezindexvalue是与Feature的几何形状相交的所有单元格的值。Featurezindex决定一个单元格在多个几何图形相交时的值。zindex 越高,优先级越高。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
geom1 = box(0, 0, 15, 15)
geom2 = box(100, 26, 109, 208)
geom3 = box(610, 215, 1000, 500)

cell_value1 = gps.CellValue(value=1, zindex=1)
cell_value2 = gps.CellValue(value=2, zindex=2)
cell_value3 = gps.CellValue(value=3, zindex=3)

# Will not be selected if feature2 and/or feature3 also intersects target cell
feature1 = gps.Feature(geometry=geom1, properties=cell_value1)

# Will not be selected if feature3 also intersects target cell
feature2 = gps.Feature(geometry=geom2, properties=cell_value2)

# Will always be selected
feature3 = gps.Feature(geometry=geom3, properties=cell_value3)