微软行星云计算Planetary Computer——重投影和重采样（方法1）

如何将云优化GeoTIFFs（COG）重新投影到不同的坐标参考系（CRS）？重投数据通常还包括将光栅数据重新采样到不同的分辨率，这需要将多个像素值组合到一起（当下采样--或重新采样到较低的分辨率/较大的单元尺寸时）或从现有的像素中插值新的像素（当上采样--或重新采样到较高的分辨率/较小的单元尺寸时）。你还将学习如何使用不同的重采样算法来调整光栅的像素大小。

本教程涵盖了几种再现的方法。选择哪种方法主要取决于两个因素。

当你加载数据时，你是否已经知道你的目标是什么CRS

你是否有用于你所处理的数据的STAC项目？

一般来说，在加载数据的同时进行重投影和重采样会更有效率，行星计算机数据目录中的所有遥感数据集都包括STAC项目。

本教程包括以下四种对地理数据进行重投影和重采样的方法。

方法1：在使用stackstac加载时对大量的COG进行重投和重采样。这种方法需要STAC项目。

方法2：在使用rioxarray和rasterio加载时，使用WarpedVRT对单个COG进行重投影和重采样。这种方法在没有STAC项目的情况下也可以使用，但是不能扩展到大量的COG。

方法3：用stackstac.reproject_array对内存中（可能是分布式的）DataArray进行重投和重新采样。这个选项需要更多的内存，一般来说效率不如选项1，但你在加载数据时不需要知道目标CRS。

方法4: 使用rioxarray和rasterio对内存中（单机）的DataArray进行重投和重新采样。这个选项只适用于单机，并且比选项2需要更多的内存，但是在加载数据时，你不需要知道目标CRS。

这些重投方法是按照一般推荐的顺序排列的。然而，某些方法可能无法使用，这取决于你在加载数据时是否知道你的目标CRS以及你是否有STAC项目。

在开始方法前，首先要进行库的导入：

import planetary_computer as pc
import pystac_client
import numpy as np
import xarray as xr
import stackstac
import affine
import pyproj
import rioxarray
import rasterio
from rioxarray.rioxarray import _make_coords
from rasterio.vrt import WarpedVRT
import xrspatial.multispectral as ms
import matplotlib.pyplot as plt

准备工作。创建一个本地Dask集群

在本教程中，将使用一个小型数据集。创建一个本地Dask集群，使用你机器的所有核心来并行处理数据。

from dask.distributed import Client
client = Client()
print(f"/proxy/{client.scheduler_info()['services']['dashboard']}/status")

/proxy/8787/status

为了跟踪计算的进展，你可以通过前一个单元的输出的URL访问Dask仪表板。

本教程所关注的区域位于美国科罗拉多州，以Cochetopa Dome和Sawtooth Mountain为中心。将使用Sentinel-2 Level-2A数据。关于通过STAC API访问数据的更多信息，请参见从STAC API中读取数据。

stac = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1/")
search = stac.search(
    collections=["sentinel-2-l2a"],
    intersects={"type": "Point", "coordinates": [-106.714900, 38.22679]},
    query={"eo:cloud_cover": {"lt": 1}},
    datetime="2019-09-01/2019-09-30",
)
item = next(search.get_items())
signed_item = pc.sign(item)

方法1： 使用stackstac

选定的场景被存储在EPSG:32613坐标系中。让我们通过向stackstac.stack函数提供目标epsg编号，将其重新投射到一个新的CRS中。要调整数据的分辨率，还需要向函数的可选分辨率和重采样参数提供数值。

通过 stackstac.stack，你可以使用任何由 EPSG 代码定义的 CRS。在这个例子中，你使用stackstac.stack函数来重新投射到Lambert Cylindrical投影（epsg:6933）。兰伯特圆柱投影是一种保留面积测量的等面积投影。在地球表面具有相同大小的区域在地图上具有相同的大小。然而，形状、角度和比例尺一般会被扭曲。

这个例子使用了双线性重采样方法。重采样也发生在stackstac.stack函数中。由于这个函数内部使用了rasterio，你可以使用rasterio.enums.Resampling类中的任何算法。要使用一个特定的算法，把它作为重采样参数的一个值来传递。stackstac.stack中默认的重采样方法是rasterio.enums.Resampling.nearest。

scene_data = (
    stackstac.stack(
        [signed_item.to_dict()],
        epsg=6933,
        resampling=rasterio.enums.Resampling.bilinear,
        resolution=100,  # resolution in the output CRS’s units
        assets=["B04", "B03", "B02"],  # red, green, blue bands
        chunksize=2048,
    )
    .isel(time=0)
    .persist()
)
scene_data

<xarray.DataArray 'stackstac-70964429577d6f8e4004178ae459a9c1' (band: 3, y: 995, x: 1235)>

dask.array<getitem, shape=(3, 995, 1235), dtype=float64, chunksize=(1, 995, 1235), chunktype=numpy.ndarray>

Coordinates: (12/46)

   time                                     datetime64[ns] 2019-09-24T17:50:...

   id                                       <U54 'S2B_MSIL2A_20190924T175049...

 * band                                     (band) <U3 'B04' 'B03' 'B02'

 * x                                        (x) float64 -1.035e+07 ... -1.02...

 * y                                        (y) float64 4.591e+06 ... 4.491e+06

   s2:product_type                          <U7 'S2MSI2A'

   ...                                       ...

   proj:transform                           object {0.0, 300000.0, 10.0, 430...

   proj:bbox                                object {300000.0, 4190220.0, 430...

   common_name                              (band) <U5 'red' 'green' 'blue'

   center_wavelength                        (band) float64 0.665 0.56 0.49

   full_width_half_max                      (band) float64 0.038 0.045 0.098

   epsg                                     int64 6933

Attributes:

   spec:        RasterSpec(epsg=6933, bounds=(-10353400, 4491300, -10229900,...

   crs:         epsg:6933

   transform:   | 100.00, 0.00,-10353400.00|\n| 0.00,-100.00, 4590800.00|\n|...

   resolution:  100

xarray.DataArray

'stackstac-70964429577d6f8e4004178ae459a9c1'

• band: 3
• y: 995
• x: 1235
• dask.array<chunksize=(1, 995, 1235), meta=np.ndarray>

• 
  
• Coordinates: (46)

• time
  ()
  datetime64[ns]
  2019-09-24T17:50:49.024000

array('2019-09-24T17:50:49.024000000', dtype='datetime64[ns]')

• id
  ()
  <U54
  'S2B_MSIL2A_20190924T175049_R141...

array('S2B_MSIL2A_20190924T175049_R141_T13SCC_20201004T164728',

     dtype='<U54')

• band
  (band)
  <U3
  'B04' 'B03' 'B02'

array(['B04', 'B03', 'B02'], dtype='<U3')

• x
  (x)
  float64
  -1.035e+07 ... -1.023e+07

array([-10353400., -10353300., -10353200., ..., -10230200., -10230100.,

      -10230000.])

• y
  (y)
  float64
  4.591e+06 4.591e+06 ... 4.491e+06

array([4590800., 4590700., 4590600., ..., 4491600., 4491500., 4491400.])

• s2:product_type
  ()
  <U7
  'S2MSI2A'

array('S2MSI2A', dtype='<U7')

• sat:orbit_state
  ()
  <U10
  'descending'

array('descending', dtype='<U10')

• sat:relative_orbit
  ()
  int64
  141

array(141)

• s2:processing_baseline
  ()
  <U5
  '02.12'

array('02.12', dtype='<U5')

• eo:cloud_cover
  ()
  float64
  0.04679

array(0.046785)

• platform
  ()
  <U11
  'Sentinel-2B'

array('Sentinel-2B', dtype='<U11')

• proj:epsg
  ()
  int64
  32613

array(32613)

• ........
• 
  
• Attributes: (4)spec :
  RasterSpec(epsg=6933, bounds=(-10353400, 4491300, -10229900, 4590800), resolutions_xy=(100, 100))
  crs :
  epsg:6933
  transform :
  | 100.00, 0.00,-10353400.00| | 0.00,-100.00, 4590800.00| | 0.00, 0.00, 1.00|
  resolution :
  100
  
• 处理结果我就不逐一展示了

为了生成重投影和重采样数据的预览，用xrspatil.multispectral.true_color函数创建一个真彩色图像。

提供以下参数。

你的场景数据中的红、绿、蓝三个波段

用于图像名称的可选字符串名称

注意：要调整对比度和亮度，请提供两个可选的数字参数c（对比度）和th（亮度）的值。本教程中的所有图像预览都使用该函数的默认值c=10和th=0.125。

结果展示：

# visualize selected scene
cylindrical_img = ms.true_color(*scene_data, name="epsg=6933")
cylindrical_img.plot.imshow(figsize=(8, 8));

下次将介绍改变投影和重采样的方法2