搜 索

从多幅时序栅格构建时空立方体的两种方案

  • 897阅读
  • 2022年04月19日
  • 0评论
首页 / 学习区 / 正文

本文介绍了两种从多幅时序栅格构建时空立方体的方法。这适用对长时序的栅格场数据进行时空模式挖掘的准备工作,以在 ArcGIS Pro 里进行时空模式挖掘

使用xarray构建多维栅格

环境准备

  • rasterio 是一个高性能的、方便易用的、专门读写栅格文件的 Python 库;
  • dask 与 numpy Array 与 pandas DataFrame 有机集成、支持在大于内存的数据集上提供多核、并行执行;
  • xarray 可以处理 N 维数组和数据集,能够读写 NetCDF、HDF 和 GRIB 等格式的数据

由于 rasterio 库有些难以安装,建议在 conda 中新开一个环境,而后安装各个包。

conda create --name nc
conda activate nc
conda install rasterio dask xarray

将多幅栅格直接合并至netcdf

这里借鉴了Stackoverflow上的回答,做了一些小的补充。实际上可以使用 xr.stack()或其它方法。

import pandas as pd
import glob
from osgeo import gdal
import xarray as xr

## 1、构建时间索引
def get_time_idx_from_tifs(tifs_path):
    '''Create a DatetimeIndex with tifs naming like `20150520.tif` '''
    return pd.DatetimeIndex([pd.Timestamp(i[-12:-4]) for i in tifs_path])

file_list = sorted(glob.glob(r"D:\Desktop\data\*.tif"))
t_idx = get_time_idx_from_tifs(file_list)
t_var = xr.Variable('time', t_idx)

## 2、将多个 tif 文件按顺序读入一个 dask Araay
ds = gdal.Open(file_list[0])
x, y = ds.RasterXSize, ds.RasterYSize
chunks = {'x': x, 'y': y, 'band': 1}
da = xr.concat([xr.open_rasterio(i, chunks=chunks) for i in file_list], dim=t_var)

## 3、将 DataArray 转换为 Dataset,以重命名变量名
ds = da.to_dataset(name='New_area')

## 4、导出 Dataset 至 netcdf
ds.to_netcdf("result.nc")

原回答的解决方案里,导出的变量会是默认的__xarray_dataarray_variable__,很不美观,原因可以参考这个链接。可以先将DataArray转换为Dataset,这样就可以修改导出 nc 文件的变量名了。

不足之处

经实践,如此生成的 .nc 文件可以在 ArcGIS Pro 中转换为 .crf 格式的多维栅格,也可以在其中进行一系列分析,但是无法再进一步生成 .nc 格式的时空立方体。在我的实践中,ArcGIS Pro 会在生成时空立方体时会报错 ERROR 110295,说我“输入的多维栅格图层不包含有效的时间维度”。但我的多维栅格可以正确滑动时间滑块,也可以正确生成时态图,检查、检索许久仍未能解决这个问题。

所以如果有遇到同样问题的 uu 解决了这个问题,可以告诉我一下该如何解决咩?

通过多维镶嵌数据集构建多维栅格

经实践,上述方案需要稍微配置下环境,结果也可能会报错,而下面这个方法不需要安装任何其它的包,只需要组合使用 ArcGIS Pro 里的若干工具。通过多维镶嵌数据集构建多维栅格的方案如下:

1、Create Mosaic Dataset

这里以土地利用数据为例,在创建时指定了像素类型为 8-bit unsigned (0-255)。需要注意指定好投影坐标系。

image-20220414200258748

2、Add Rasters To Mosaic Dataset

上一步创建了一个空的镶嵌数据集,接下来需要将多个年份的 tif 文件按顺序填入其中,注意统一投影。此外,建议将各栅格文件的名称就设置为年份日期,便于后续处理!

image-20220414200601020

3、Convert Time Field

将多幅栅格加入镶嵌数据集后,数据集下的的 Footprint 属性表中会加载一系列信息,需要将文件名称字段转换为 Date 格式,以便于多维栅格读取。

image-20220414202615397

4、Build Multidimensional Info

为镶嵌数据集写入多维信息!注意在维度 Dimension 一栏写时间信息,而变量 Variable 一栏则写所分析的变量名(如温度、人口数等)。我这里出于演示仅用了属性表中的 Dataset 字段,这里可以自行调整!

image-20220414203758247

5、Copy Raster

这一步很重要,将多维镶嵌数据集 Mosaic 转为多维云栅格 .crf,这里也支持导出.nc.zarr等格式。

image-20220414204412597

复制完成后,可以右击属性表中确认多维栅格的相关信息。

image-20220414203645959

6、Create Space Time Cube From Multidimensional Raster Layer

最后就可以从多维栅格创建时空立方体了,大功告成!

image-20220414204959376

小结

经实践,构建时空立方体的稳定方案是:多幅时序栅格→多维镶嵌数据集→多维栅格→时空立方体。

评论区
暂无评论
avatar