GeoPySpark处理各类栅格数据

GeoPySpark处理各类栅格数据

在开始之前,我们首先下载样例数据和导入模块包:

1
curl -o /Users/xujavy/Documents/Work/data/jupyter_data/cropped.tif https://s3.amazonaws.com/geopyspark-test/example-files/cropped.tif

1
2
3
# geopyspark需要设置SPARK_HOME
import os
os.environ["SPARK_HOME"] = '/Users/xujavy/Documents/Work/geoserver/installpackages/spark-2.2.3-bin-hadoop2.7'
1
2
3
4
5
6
7
8
9
10
import datetime
import numpy as np
import pyproj
import geopyspark as gps

from pyspark import SparkContext
from shapely.geometry import box, Point

conf = gps.geopyspark_conf(master="local[*]", appName="layers")
pysc = SparkContext(conf=conf)

如何在GeoPySpark中存储和表示数据?

GeoPySpark中使用的所有数据都存储在RDD中。因此,理解GeoPySpark如何在整个库中存储、表示和使用这些RDD非常重要。

GeoPySpark不使用PySpark RDDs,而是使用Python类作为Scala类的包装器,这些类包含并使用Scala RDD。具体来说,这些包装器类是RasterLayerTiledRasterLayer层,后面将更详细地讨论它们。

图层不仅仅是RDD

我们将Python包装器类称为层而不是RDDs,原因有两个:首先,RasterLayerTiledRasterLayer实际上都没有扩展PySpark的RDD类;但更重要的是,这些类包含比RDD更多的信息。当我们提到“层”时,我们指的是RDD及其属性。

GeoPySpark层包含的RDDs包含类型为(K, V)的元组,其中K表示键,V表示值。V始终是Tile,但是K的不同取决于包装器类和数据本身的性质。更多信息请见下文。

RasterLayer

RasterLayer类处理未处理的数据,也就是说,该图层的元素没有被规范化为一个统一的布局。每个栅格元素可能有不同的分辨率或尺寸;组成栅格的范围不需要遵循任何有序的模式。本质上,RasterLayer存储”raw”数据,其主要目的是充当获取遵循指定布局的瓦片数据的路径上的中转站。

RasterLayer对象包含RDDs分别为ProjectedExtentTemporalProjectedExtent的键类型K,该层类型为SPATIALPACETIME

TiledRasterLayer

TiledRasterLayerRasterLayer的补充,用于存储瓦片数据。瓦片数据已符合一定的布局,这意味着它已定期取样,并已被分割成大小均匀、不重叠的片段,可以合理地建立索引。让数据处于这种状态的好处是,很容易使用它。例如,通过这个类,用户将能够执行映射代数、创建金字塔和保存层。下面是这些操作的定义和具体示例。

对于TiledRasterLayer, K是SpatialKeySpaceTimeKey

RasterLayer操作

创建RasterLayers

有两种方法可以创建RasterLayer:(1)通过本地文件系统、S3或HDFS读取geotiff;(2)来自现有的PySpark RDD。

来自PySpark RDDs

从PySpark RDD创建RasterLayer第一个参数是from_numpy_rdd() 类方法。这一步可能有点复杂,因为它要求以特定的方式格式化PySpark RDD中的数据。

下面的示例从元组构造RDD。第一个参数是ProjectedExtent,因为我们决定将数据空间化。如果我们处理的是时空数据,那么TemporalProjectedExtent将是第一个参数。Tile总是元组的第二个元素。

1
2
3
4
5
6
7
8
9
arr = np.ones((1, 16, 16), dtype='int')
tile = gps.Tile.from_numpy_array(numpy_array=np.array(arr), no_data_value=-500)

extent = gps.Extent(0.0, 1.0, 2.0, 3.0)
projected_extent = gps.ProjectedExtent(extent=extent, epsg=3857)

rdd = pysc.parallelize([(projected_extent, tile), (projected_extent, tile)])
multiband_raster_layer = gps.RasterLayer.from_numpy_rdd(layer_type=gps.LayerType.SPATIAL, numpy_rdd=rdd)
multiband_raster_layer
来自GeoTiffs

geopyspark.geotrellis中的get() 函数,可以从geotiff创建一个RasterLayer实例。这些文件可以位于本地文件系统、HDFS或S3上。在本例中,本地读取带有空间数据的GeoTiff。

1
2
raster_layer = gps.geotiff.get(layer_type=gps.LayerType.SPATIAL, uri="file:///Users/xujavy/Documents/Work/data/jupyter_data/cropped.tif")
raster_layer

使用RasterLayer

下一节将介绍RasterLayer方法。应该注意的是,并不是该类中包含的所有方法都将被覆盖。在GeoPySpark中的可视化数据中可以找到更多关于处理层内容可视化的方法的信息。

转换成一个Python RDD

通过使用to_numpy_rdd()RasterLayer将被序列化为Python RDD。这将把每个元组中的第一个值转换为ProjectedExtentTemporalProjectedExtent,第二个值转换为Tile

1
2
python_rdd = raster_layer.to_numpy_rdd()
python_rdd
1
python_rdd.first()
SpaceTime Layer 转换为 Spatial Layer

如果使用的是spatial-temporal图层,并且希望将其转换为空间层,那么可以使用to_spatial_layer() 方法。这将通过将TemporalProjectedExtent转换为ProjectedExtent来更改层中RDD的键。

1
2
3
4
5
6
7
8
9
# Creating the space time layer
instant = datetime.datetime.now()
temporal_projected_extent = gps.TemporalProjectedExtent(extent=projected_extent.extent,
epsg=projected_extent.epsg,
instant=instant)

space_time_rdd = pysc.parallelize([temporal_projected_extent, tile])
space_time_layer = gps.RasterLayer.from_numpy_rdd(layer_type=gps.LayerType.SPACETIME, numpy_rdd=space_time_rdd)
space_time_layer
1
2
# Converting the SpaceTime layer to a Spatial layer
space_time_layer.to_spatial_layer()
收集元数据

图层的Metadata包含图层中值的信息。这些数据包含该图层的的布局、投影和范围。

collect_metadata()将返回符合给定layout图层的Metadata

1
2
3
# Collecting Metadata with the default LocalLayout()
metadata = raster_layer.collect_metadata()
metadata
1
2
# Collecting Metadata with the default GlobalLayout()
raster_layer.collect_metadata(layout=gps.GlobalLayout())
1
2
3
4
5
6
# Collecting Metadata with a LayoutDefinition
extent = gps.Extent(0.0, 0.0, 33.0, 33.0)
tile_layout = gps.TileLayout(2, 2, 256, 256)
layout_definition = gps.LayoutDefinition(extent, tile_layout)

raster_layer.collect_metadata(layout=layout_definition)
重定义投影

reproject() 将把图层内栅格数据的原始投影更改为给定的target_crs投影。此方法不会对超过瓦片边界的区域进行采样。

1
2
# The CRS of the layer before reprojecting
metadata.crs
1
2
# The CRS of the layer after reprojecting
raster_layer.reproject(target_crs=3857).collect_metadata().crs
将Tile数据转换到一个Layout

tile_to_layout()RasterLayer的栅格瓦片格式化为给定的布局。这个切片的结果是一个新的TiledRasterLayer实例。此输出包含与源RasterLayer相同的数据,但是,其中包含的信息现在将根据给定的布局进行组织。

在此步骤中,还可以重新投影RasterLayer。这可以通过指定要重新映射到的target_crs来实现。使用此方法进行重新投影产生的结果与重新投影方法返回的结果不同。后者不会超出图层内栅格的边界进行采样,而前者会。这一点非常重要,因为全局布局需要超越栅格的边界进行采样。

  1. 从Metadata

    创建一个包含给定Metadata布局的TiledRasterLayer

    注意: 如果指定的target_crs与元数据中的不同,那么将抛出一个错误。

    1
    raster_layer.tile_to_layout(layout=metadata)
  2. 从LayoutDefinition

    1
    raster_layer.tile_to_layout(layout=layout_definition)
  3. 从LocalLayout

    1
    raster_layer.tile_to_layout(gps.LocalLayout())
  4. 从GlobalLayout

    1
    2
    tiled_raster_layer = raster_layer.tile_to_layout(gps.GlobalLayout())
    tiled_raster_layer
  5. 从一个TiledRasterLayer
    可以将RasterLayertile的tile与TiledRasterLayout相同的布局。

    注意:如果指定的target_crs与其他层的不同,那么将抛出一个错误。

    1
    raster_layer.tile_to_layout(layout=tiled_raster_layer)

TiledRasterLayer

创建TiledRasterLayers

我们将介绍一个TiledRasterLayer的初始化方法from_numpy_rdd。然而,还有其他方法可以创建这个类。这些额外的创建策略可以在[map algebra guide]中找到。

从PySpark RDD

RasterLayer一样,可以使用from_numpy_rdd()RDDs创建TiledRasterLayer。但是,不同的是,元数据还必须在初始化期间传入。这使得以这种方式创建TiledRasterLayer变得更加困难。

下面的示例从元组构造RDD。第一个元素是SpatialKey,因为我们决定将数据空间化。

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
data = np.zeros((1, 512, 512), dtype='float32')
tile = gps.Tile.from_numpy_array(numpy_array=data, no_data_value=-1.0)
instant = datetime.datetime.now()

layer = [(gps.SpaceTimeKey(row=0, col=0, instant=instant), tile),
(gps.SpaceTimeKey(row=1, col=0, instant=instant), tile),
(gps.SpaceTimeKey(row=0, col=1, instant=instant), tile),
(gps.SpaceTimeKey(row=1, col=1, instant=instant), tile)]

rdd = pysc.parallelize(layer)

extent = gps.Extent(0.0, 0.0, 33.0, 33.0)
layout = gps.TileLayout(2, 2, 512, 512)
bounds = gps.Bounds(gps.SpaceTimeKey(col=0, row=0, instant=instant), gps.SpaceTimeKey(col=1, row=1, instant=instant))
layout_definition = gps.LayoutDefinition(extent, layout)

metadata = gps.Metadata(
bounds=bounds,
crs='+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 ',
cell_type='float32ud-1.0',
extent=extent,
layout_definition=layout_definition)

space_time_tiled_layer = gps.TiledRasterLayer.from_numpy_rdd(layer_type=gps.LayerType.SPACETIME,
numpy_rdd=rdd, metadata=metadata)
space_time_tiled_layer

使用TiledRasterLayers

本节将介绍在TiledRasterLayer中找到的方法。就像使用RasterLayer一样,但是也不会介绍该类中的所有方法。在GeoPySpark可视化数据中可以找到更多关于处理层内容可视化的方法的信息。

转换为一个Python RDD

通过使用to_numpy_rdd() ,基本的TiledRasterLayer将被序列化为Python RDD。这将把每个元组中的所有第一个值转换为SpatialKeySpaceTimeKey,第二个值转换为Tile

1
2
python_rdd = tiled_raster_layer.to_numpy_rdd()
python_rdd.first()
SpaceTime Layer转换为Spatial Layer

如果使用的是一个时空层,并且希望将其转换为一个空间层,那么可以使用to_spatial_layer() 方法。这将通过将SpaceTimeKey转换为SpatialKey来更改层中RDD的键。

1
2
# Converting the SpaceTime layer to a Spatial layer
space_time_tiled_layer.to_spatial_layer()
  1. 重新分区

    虽然不是RDD,但TiledRasterLayer确实包含底层RDD,因此可以使用repartition() 方法对其进行重新分区。

    1
    2
    # Repartition the internal RDD to have 120 partitions
    tiled_raster_layer.repartition(num_partitions=120)
  2. 查找

    如果该层中存在感兴趣的特定平铺,可以使用lookup() 方法将其作为Tile检索。

    1
    2
    3
    4
    min_key = tiled_raster_layer.layer_metadata.bounds.minKey

    # Retrieve the Tile that is located at the smallest column and row of the layer
    tiled_raster_layer.lookup(col=min_key.col, row=min_key.row)
  3. 掩码

    通过使用mask() 方法,可以使用一个或多个形状几何图形对·TiledRasterRDD·进行掩码。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    layer_extent = tiled_raster_layer.layer_metadata.extent

    # Polygon to mask a region of the layer
    mask = box(layer_extent.xmin,
    layer_extent.ymin,
    layer_extent.xmin + 20,
    layer_extent.ymin + 20)

    tiled_raster_layer.mask(geometries=mask)
    1
    2
    3
    4
    5
    6
    7
    mask_2 = box(layer_extent.xmin + 50,
    layer_extent.ymin + 50,
    layer_extent.xmax - 20,
    layer_extent.ymax - 20)

    # Multiple Polygons can be given to mask the layer
    tiled_raster_layer.mask(geometries=[mask, mask_2])
  4. 正常化

    normalize() 将对图层中的数据进行线性转换,使所有值都落在给定的范围内。

    1
    2
    # Normalizes the layer so that the new min value is 0 and the new max value is 60000
    tiled_raster_layer.normalize(new_min=0, new_max=60000)
  5. 创建金字塔

    在为TMS服务器使用图层时,重要的为图层创建金字塔。也就是说,我们创建一个覆盖相同地理范围的细节层次结构,而金字塔的每一层使用的像素是下一层的四分之一。这使我们可以放大和缩小时,图层显示没有使用无关的细节。pyramid() 方法将生成一个金字塔实例,该实例将在其中包含多个TiledRasterLayer。每个层对应一个缩放级别,级别的数量取决于源层的zoom_level金字塔的最大级别为源层的zoom_level,最小级别为0。

    1
    2
    # This creates a Pyramid with zoom levels that go from 0 to 11 for a total of 12.
    tiled_raster_layer.pyramid()
  6. 重投影

    这类似于RasterLayerreproject方法,在该方法中,重新投影将不会采样超过瓦片的边界。这意味着tile的布局将会改变,因此它们将呈现LocalLayout而不是GlobalLayout。因此,无论TiledRasterLayerzoom_level是什么,都将更改为0,因为所表示的区域只更改为tile。

    1
    2
    # The zoom_level and crs of the TiledRasterLayer before reprojecting
    tiled_raster_layer.zoom_level, tiled_raster_layer.layer_metadata.crs
    1
    2
    3
    4
    reprojected_tiled_raster_layer = tiled_raster_layer.reproject(target_crs=3857)

    # The zoom_level and crs of the TiledRasterLayer after reprojecting
    reprojected_tiled_raster_layer.zoom_level, reprojected_tiled_raster_layer.layer_metadata.crs
  7. 缝合

    使用stitch()TiledRasterLayer内的所有Tile缝合在一起,得到一个Tile。这只能在空间层中完成,如果层中包含的数据很大,则不建议这样做,因为它可能会由于结果块的大小而导致崩溃。

    1
    2
    # Creates a Tile with an underlying numpy array with a size of (1, 6144, 1536).
    tiled_raster_layer.stitch().cells.shape
  8. 保存一个缝合的图层

    save_stitched() 方法既可以缝合,也可以将层保存为GeoTiff。

    1
    2
    # Saves the stitched layer to /tmp/stitched.tif
    tiled_raster_layer.save_stitched(path='./stitched.tif')

    也可以指定缝制层时要保存的区域。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    layer_extent = tiled_raster_layer.layer_metadata.layout_definition.extent

    # Only a portion of the stitched layer needs to be saved, so we will create a sub Extent to crop to.
    sub_exent = gps.Extent(xmin=layer_extent.xmin + 10,
    ymin=layer_extent.ymin + 10,
    xmax=layer_extent.xmax - 10,
    ymax=layer_extent.ymax - 10)

    tiled_raster_layer.save_stitched(path='./cropped-stitched.tif', crop_bounds=sub_exent)
    1
    2
    3
    4
    # In addition to the sub Extent, one can also choose how many cols and rows will be in the saved in the GeoTiff.
    tiled_raster_layer.save_stitched(path='./cropped-stitched-2.tif',
    crop_bounds=sub_exent,
    crop_dimensions=(1000, 1000))
  9. 将tile数据转换为Layout

    这类似于RasterLayertile_to_layout方法,除了一个重要的细节。如果在包含zoom_levelTiledRasterLayer上执行tile_to_layout() ,该zoom_level可能会丢失或更改,这取决于选择的layouttarget_crs。因此,在重新排列TiledRasterLayer时,一定要记住这一点。

    1
    2
    # Original zoom_level of the source TiledRasterLayer
    tiled_raster_layer.zoom_level
    1
    2
    # zoom_level will be lost in the resulting TiledRasterlayer
    tiled_raster_layer.tile_to_layout(layout=gps.LocalLayout())
    1
    2
    # zoom_level will be changed in the resulting TiledRasterLayer
    tiled_raster_layer.tile_to_layout(layout=gps.GlobalLayout(), target_crs=3857)
    1
    2
    # zoom_level will reamin the same in the resulting TiledRasterLayer
    tiled_raster_layer.tile_to_layout(layout=gps.GlobalLayout(zoom=11))
  10. 获取点值

    get_point_values()接受shapely.geometry.Point。返回图层中给定点的值。返回的值的数量取决于值的波段,因为每个波段将有一个值。

    也可以将ResampleMethod传递给此方法,但不支持所有方法。下面是所有可以用来计算点值的ResampleMethods:

    • ResampleMethod.NEAREST_NEIGHBOR
    • ResampleMethod.BILINEAR
    • ResampleMethod.CUBIC_CONVOLUTION
    • ResampleMethod.CUBIC_SPLINE

    从空间图层获取点值

    当在LayerTypeSPATIAL的图层上使用get_point_values时,结果将成对为(shapely.geometry.Point, [float])。每个给定的Points与它所相交的值配对。

    1
    2
    3
    4
    5
    # Creating the points
    extent = tiled_raster_layer.layer_metadata.extent

    p1 = Point(extent.xmin, extent.ymin + 0.5)
    p2 = Point(extent.xmax , extent.ymax - 1.0)

    从[shapely.geometry.Point]获取点值

    points是一个[shapely.geometry.Point], 将返回一个[(shapely.geometry.Point, [float])].

    1
    tiled_raster_layer.get_point_values(points=[p1, p2])

    从{k: shapely.geometry.Point}获取点值

    points是一个{k: shapely.geometry.Point}, 将返回一个{k: (shapely.geometry.Point, [float])}

    1
    tiled_raster_layer.get_point_values(points={'point 1': p1, 'point 2': p2})

    从时空图层获取点值

    当在LayerTypeSPACETIME的图层上使用get_point_values时,结果将成对为(shapely.geometry.Point, [(datetime.datetime, [float])])。每个给定的Points将与一个元组列表配对,其中包含与之相交的值以及这些值的相应时间戳。

    1
    2
    3
    4
    st_extent = space_time_tiled_layer.layer_metadata.extent

    p1 = Point(st_extent.xmin, st_extent.ymin + 0.5)
    p2 = Point(st_extent.xmax , st_extent.ymax - 1.0)

    从[shapely.geometry.Point]获取点值

    points是一个[shapely.geometry.Point], 将返回一个[(shapely.geometry.Point, [(datetime.datetime, [float])])].

    1
    space_time_tiled_layer.get_point_values(points=[p1, p2])

    从{k: shapely.geometry.Point}获取点值

    points是一个{k: shapely.geometry.Point}, 将返回一个[(shapely.geometry.Point, [(datetime.datetime, [float])])].

    1
    space_time_tiled_layer.get_point_values(points={'point 1': p1, 'point 2': p2})
  11. 聚合每个单元格的值
    aggregate_by_cell() 将为每个单元计算每个键的所有值的聚合汇总。因此,如果图层中有相同键的多个副本,那么结果图层将只包含该键的一个实例,其对应的值是共享该键的所有值的汇总。

    不支持所有操作。下面这些可以在aggregate_by_cell中使用:

    • Operation.SUM
    • Operation.MIN
    • Operation.MAX
    • Operation.MEAN
    • Operation.VARIANCE
    • Operation.STANDARD_DEVIATION
    1
    2
    3
    4
    5
    6
    unioned_layer = gps.union(layers=[tiled_raster_layer, tiled_raster_layer + 1])
    #Sum the values of the unioned_layer
    unioned_layer.aggregate_by_cell(operation=gps.Operation.SUM)

    # Get the max value for each cell
    unioned_layer.aggregate_by_cell(operation=gps.Operation.MAX)

一般的方法

RasterLayerTiledRasterLayer中都有一些方法。这些方法倾向于执行更一般的分析/任务,因此它们适合这两个类。下一节将介绍这些方法。

注意:在下面的示例中,将同时使用RasterLayerTiledRasterLayer。然而,它们很容易被其他类替换。

结合图层

要将多层的内容组合在一起,可以使用union()方法。这将生成一个新的包含来自给定层的所有元素的RasterLayerTiledRasterLayer

注意: 生成的层可以包含重复的键。

1
gps.union(layers=[tiled_raster_layer, tiled_raster_layer])

选择波段的一个子波段

要选择要使用的特定波段,bands()方法将采用单个波段索引或波段索引集合,并将该子集作为新的RasterLayerTiledRasterLayer返回。

注意: 如果在一个大数据集的两个子带之间执行操作,可能会有很高的性能成本。因此,如果您正在处理大量数据,那么建议在读取它们之前进行波段选择。

1
2
# Selecting the second band from the layer
multiband_raster_layer.bands(1)

1
2
# Selecting the first and second bands from the layer
multiband_raster_layer.bands([0, 1])

两个或更多层的波段组合

combine_bands() 方法将连接两个或多个层之间共享键值的波段。因此,结果图层将为每个共享键值包含一个新的Tile, Tile将包含来自给定层的所有波段。

图层传递到combine_bands的顺序很重要。其中,结果值的波段将根据其各自层的位置进行排序。

1
2
3
4
5
6
# Setting up example RDD
twos = np.ones((1, 16, 16), dtype='int') + 1
twos_tile = gps.Tile.from_numpy_array(numpy_array=np.array(twos), no_data_value=-500)

twos_rdd = pysc.parallelize([(projected_extent, twos_tile)])
twos_raster_layer = gps.RasterLayer.from_numpy_rdd(layer_type=gps.LayerType.SPATIAL, numpy_rdd=twos_rdd)

1
2
3
# The resulting values of the layer will have 2 bands: the first will be all ones,
# and the last band will be all twos
gps.combine_bands(layers=[multiband_raster_layer, twos_raster_layer])
1
2
3
# The resulting values of the layer will have 2 bands: the first will be all twos and the
# other band will be all ones
gps.combine_bands(layers=[twos_raster_layer, multiband_raster_layer])

收集图层的键值

要收集一个层的所有键,使用collect_keys方法。

1
2
3
4
5
6
7
8
# Returns a list of ProjectedExtents
multiband_raster_layer.collect_keys()

# Returns a list of a SpatialKeys
tiled_raster_layer.collect_keys()

# Returns a list of SpaceTimeKeys
space_time_tiled_layer.collect_keys()

按时间过滤图层

使用filter_by_times方法将生成一个值位于给定时间间隔内的层。

按精确时间过滤

一个datetime.datetime实例可用于过滤该层。如果是这样,那么只保留与给定时间的精确匹配。

1
space_time_layer.filter_by_times(time_intervals=[instant])

按时间间隔过滤

也可以给出不同的时间间隔,任何时间落在时间跨度内的键都将保存在层中。

1
2
3
4
5
6
7
8
9
10
end_date_1 = instant + datetime.timedelta(days=3)
end_date_2 = instant + datetime.timedelta(days=5)

# Will filter out any value whose key does not fall in the range of
# instant and end_date_1
space_time_layer.filter_by_times(time_intervals=[instant, end_date_1])

# Will filter out any value whose key does not fall in the range of
# instant and end_date_1 OR whose key does not match end_date_2
space_time_layer.filter_by_times(time_intervals=[instant, end_date_1, end_date_2])

转换光栅单元格的数据类型

convert_data_type 方法将图层的栅格中的单元格类型转换为新的数据类型。在此转换期间还可以设置noData值,如果没有设置,则生成的栅格数据将没有noData值。

1
2
# The data type of the cells before converting
metadata.cell_type

1
2
# Changing the cell type to int8 with a noData value of -100.
raster_layer.convert_data_type(new_type=gps.CellType.INT8, no_data_value=-100).collect_metadata().cell_type
1
2
# Changing the cell type to int32 with no noData value.
raster_layer.convert_data_type(new_type=gps.CellType.INT32).collect_metadata().cell_type

重新分类单元的值

根据给定的value_mapclassification_strategy对单元格值进行reclassify。除了这两个参数外,还需要给出单元格的data_type, 格式是int或者float

1
2
# Values of the first tile before being reclassified
multiband_raster_layer.to_numpy_rdd().first()[1]

1
2
3
4
5
# Change all values greater than or equal to 1 to 10
reclassified = multiband_raster_layer.reclassify(value_map={1: 10},
data_type=int,
classification_strategy=gps.ClassificationStrategy.GREATER_THAN_OR_EQUAL_TO)
reclassified.to_numpy_rdd().first()[1]

将图层的值合并在一起

通过使用merge方法,将图层中共享一个Key值的所有值合并在一起,形成一个新的单值。这是通过用一个值替换另一个值的单元格来实现的。然而,并非所有的单元格(如果有的话)都可以被替换。在合并值的单元格时,将采取以下步骤来确定是否应该更改单元格的值:

  1. 如果单元格包含NoData值,则将替换该值。
  2. 如果没有设置NoData值,那么将替换vlue为0的单元格。
  3. 如果以上两个都不为真,则单元格保留其值。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# Creating the layers
no_data = np.full((1, 4, 4), -1)
zeros = np.zeros((1, 4, 4))

def create_layer(no_data_value=None):
data_tile = gps.Tile.from_numpy_array(numpy_array=no_data, no_data_value=no_data_value)
zeros_tile = gps.Tile.from_numpy_array(numpy_array=zeros, no_data_value=no_data_value)

layer_rdd = pysc.parallelize([(projected_extent, data_tile), (projected_extent, zeros_tile)])
return gps.RasterLayer.from_numpy_rdd(layer_type=gps.LayerType.SPATIAL, numpy_rdd=layer_rdd)

# Resulting layer has a no_data_value of -1
no_data_layer = create_layer(-1)

# Resutling layer has no no_data_value
no_no_data_layer = create_layer()
1
2
3
4
5
# The resulting merged value will be all zeros since -1 is the noData value
no_data_layer.merge()

# The resulting merged value will be all -1's as ``no_data_value`` was set.
no_no_data_layer.merge()

单元格映射

可以通过map_cells方法直接处理图层中的单元格。该方法接受一个函数,该函数期望一个numpy数组和一个noData值作为参数,并返回一个新的numpy数组。因此,给出的函数具有以下类型签名:
def input_function(numpy_array: np.ndarray, no_data_value=None) -> np.ndarray
然后将给定的函数应用于层中的每个Tile

注意: 为了使该方法能够运行,首先需要将内部RDD从Scala反序列化到Python,然后将其从Python序列化回Scala。因此,建议将所有函数链接在一起,以避免不必要的序列化开销。

1
2
3
4
5
def add_one(cells, _):
return cells + 1

# Mapping with a single funciton
raster_layer.map_cells(add_one)

1
2
3
4
5
def divide_two(cells, _):
return (add_one(cells) / 2)

# Chaning together two functions to be mapped
raster_layer.map_cells(divide_two)

Tile映射

map_cells一样,map_tiles将给定的函数映射到图层中的所有Tile上。它接受一个期望Tile并返回Tile的函数。因此,输入函数的类型签名为:
def input_function(tile: Tile) -> Tile
注意: 为了使该方法能够运行,首先需要将内部RDD从Scala反序列化到Python,然后将其从Python序列化回Scala。因此,建议将所有函数链接在一起,以避免不必要的序列化开销。

1
2
3
4
def minus_two(tile):
return gps.Tile.from_numpy_array(tile.cells - 2, no_data_value=tile.no_data_value)

raster_layer.map_tiles(minus_two)

计算图层的直方图

可以使用get_histogramget_class_histogram方法计算图层的直方图。这两种方法都生成直方图,但是,在生成的直方图中表示数据的方式因使用的方法不同而不同。get_histogram将生成一个值为float的直方图。而get_class_histogram返回的直方图的值为int

有关直方图类的更多信息,请参见直方图[指南]。

1
2
# Returns a Histogram whose underlying values are floats
tiled_raster_layer.get_histogram()

1
2
# Returns a Histogram whose underlying values are ints
tiled_raster_layer.get_class_histogram()

找到图层的分位数

如果希望找到没有直方图的图层的分位数断点,那么可以使用get_quantile_breaks方法。

1
tiled_raster_layer.get_quantile_breaks(num_breaks=3)

分位数为精确整数

get_quantile_breaks_exact_intget_quantile_break_breaks_exact_int 的另一个版本,它可以精确计算整数值。但是,如果层中有太多的值,那么可能会发生内存错误。

1
tiled_raster_layer.get_quantile_breaks_exact_int(num_breaks=3)

找出图层的最小值和最大值

get_min_max 方法将找到该图层的最小值和最大值。无论单元格的数据类型如何,结果总是(float, float)

1
tiled_raster_layer.get_min_max()

将图层转换为png

通过to_png_rdd方法,可以将图层中的值转换为字节形式的PNG。为了将每个值转换为PNG,需要提供一个ColorMap

除了将每个值转换为PNG外,(K, V)的结果集合将保存在Python RDD中。

1
2
3
4
5
hist = tiled_raster_layer.get_histogram()
cmap = gps.ColorMap.build(hist, 'viridis')

png_rdd = tiled_raster_layer.to_png_rdd(color_map=cmap)
png_rdd

1
2
3
4
from PIL import Image
import io
im = Image.open(io.BytesIO(png_rdd.collect()[0][1]))
im

png

将图层换为geotiff

类似于to_png_rdd,只有to_geotiff_rdd将返回Python RDD[(K, bytes)],其中字节表示一个GeoTiff。

选择一个StorageMethod

有两种不同的方式可以格式化GeoTiff:StorageMethod.STRIPEDStorageMethod.TILED。这由storage_method参数表示。默认情况下,StorageMethod.STRIPED

选择片段的大小

有两个不同的参数控制每个段的大小:rows_per_striptile_dimensions。只需要设置其中一个值,这由storage_method是什么决定。

如果storage_methodStorageMethod.STRIPEDrows_per_strip参数将要被更改。默认情况下,将计算rows_per_strip,使每个strip小于或等于8K。

如果storage_methodStorageMethod.TILEDtile_dimensions可以设置。这是一个(int, int),其中第一个值是cols的数量,第二个值是rows。默认情况下,tile_dimensions(256,256)

选择一个CompressionMethod

可以选择的两种压缩类型是:Compression.NO_COMPRESSION 或者 Compression.DEFLATE_COMPRESSION。默认情况下,压缩参数设置为Compression.NO_COMPRESSION

选择一个色彩

color_space参数确定在每个GeoTiff中应该如何组织颜色。默认情况下,它是ColorSpace.BLACK_IS_ZERO

传入一个颜色映射

可以传入ColorMap实例,以使生成的geotiff处于不同的渐变中。默认情况下,color_mapNone

1
2
3
4
5
# Creates an RDD[(K, bytes)] with the default parameters
tiled_raster_layer.to_geotiff_rdd()

# Creates an RDD whose GeoTiffs are tiled with a size of (128, 128)
tiled_raster_layer.to_geotiff_rdd(storage_method=gps.StorageMethod.TILED, tile_dimensions=(128, 128))

RDD方法

正如在TiledRasterLayer的重分区方法一节中所提到的,TiledRasterLayer的内部RDD有很多方法。这也适用于RasterLayer

下面是一个RDD列表,其中包含两个类都支持的示例。

Cache
1
raster_layer.cache()
Persist
1
2
# If no level is given, then MEMORY_ONLY will be used
tiled_raster_layer.persist()
Unpersist
1
tiled_raster_layer.unpersist()
getNumberOfPartitions
1
raster_layer.getNumPartitions()
Count
1
raster_layer.count()
isEmpty
1
raster_layer.isEmpty()