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

举报
此星光明 发表于 2022/04/14 01:17:24 2022/04/14
【摘要】 如何将云优化GeoTIFFs(COG)重新投影到不同的坐标参考系(CRS)?重投数据通常还包括将光栅数据重新采样到不同的分辨率,这需要将多个像素值组合到一起(当下采样--或重新采样到较低的分辨率/较大的单元尺寸时)或从现有的像素中插值新的像素(当上采样--或重新采样到较高的分辨率/较小的单元尺寸时)。你还将学习如何使用不同的重采样算法来...

如何将云优化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项目。

在开始方法前,首先要进行库的导入:


  
  1. import planetary_computer as pc
  2. import pystac_client
  3. import numpy as np
  4. import xarray as xr
  5. import stackstac
  6. import affine
  7. import pyproj
  8. import rioxarray
  9. import rasterio
  10. from rioxarray.rioxarray import _make_coords
  11. from rasterio.vrt import WarpedVRT
  12. import xrspatial.multispectral as ms
  13. import matplotlib.pyplot as plt

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


  
  1. from dask.distributed import Client
  2. client = Client()
  3. 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中读取数据。


  
  1. stac = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1/")
  2. search = stac.search(
  3. collections=["sentinel-2-l2a"],
  4. intersects={"type": "Point", "coordinates": [-106.714900, 38.22679]},
  5. query={"eo:cloud_cover": {"lt": 1}},
  6. datetime="2019-09-01/2019-09-30",
  7. )
  8. item = next(search.get_items())
  9. 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。

 


  
  1. scene_data = (
  2. stackstac.stack(
  3. [signed_item.to_dict()],
  4. epsg=6933,
  5. resampling=rasterio.enums.Resampling.bilinear,
  6. resolution=100, # resolution in the output CRS’s units
  7. assets=["B04", "B03", "B02"], # red, green, blue bands
  8. chunksize=2048,
  9. )
  10. .isel(time=0)
  11. .persist()
  12. )
  13. 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>

    Array Chunk
    Bytes 28.13 MiB 9.38 MiB
    Shape (3, 995, 1235) (1, 995, 1235)
    Count 3 Tasks 3 Chunks
    Type float64 numpy.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。

结果展示:


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

 

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

文章来源: blog.csdn.net,作者:此星光明2021年博客之星云计算Top3,版权归原作者所有,如需转载,请联系作者。

原文链接:blog.csdn.net/qq_31988139/article/details/123439197

【版权声明】本文为华为云社区用户转载文章,如果您发现本社区中有涉嫌抄袭的内容,欢迎发送邮件进行举报,并提供相关证据,一经查实,本社区将立刻删除涉嫌侵权内容,举报邮箱: cloudbbs@huaweicloud.com
  • 点赞
  • 收藏
  • 关注作者

评论(0

0/1000
抱歉,系统识别当前为高风险访问,暂不支持该操作

全部回复

上滑加载中

设置昵称

在此一键设置昵称,即可参与社区互动!

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。