使用 EarthPy 堆叠和裁剪tif栅格数据

用户11172986

发布于 2024-06-20 19:17:34

600

发布于 2024-06-20 19:17:34

文章被收录于专栏：气python风雨气python风雨

使用 EarthPy 堆叠和裁剪tif栅格数据

温馨提示

本文镜像：气象分析3.9

由于可视化代码过长隐藏，可点击以下链接运行Fork查看

使用 EarthPy 堆叠和裁剪tif栅格数据🔜🔜若没有成功加载可视化图，点击运行可以查看 ps:隐藏代码在【代码已被隐藏】所在行，点击所在行，可以看到该行的最右角，会出现个三角形，点击查看即可

使用 EarthPy 堆叠和裁剪tif栅格数据

注意

下面的示例将向您展示如何使用 ''es.stack（）'' 和 EarthPy 中的 ''es.crop_image（）'' 函数。

堆叠多波段影像

一些遥感数据集与每个波段一起存储在单独的文件中。然而通常，您希望在分析中同时使用所有波段。例如您需要将所有条带放在同一个文件或“堆栈”中才能绘制颜色 RGB图像。EarthPy 有一个 ''stack（）'' 函数，可让您获取一组“.tif”文件，这些文件都位于相同的空间范围、CRS 和分辨率中并将它们一起导出为一个堆叠的“.tif”文件，或者在 Python 中使用它们直接作为堆叠的 numpy 数组。

要开始使用 EarthPy ''stack（）'' 函数，请导入所需的包并创建一个要绘制的数组。下面使用颜色条将数据绘制为连续数据使用 ''plot_bands（）'' 函数

当然，小编手头没有卫星波段数据，只好拿之前的暴雨tif数据顶顶。

安装与导入库

我们将结合多个库堆叠与裁剪tif数据

In [ ]:

!pip install earthpy -i https://pypi.mirrors.ustc.edu.cn/simple/

In [1]:

import os
from glob import glob
import matplotlib.pyplot as plt
import rasterio as rio
from rasterio.plot import plotting_extent
import geopandas as gpd
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep

注意

如果您在 Windows 系统上运行此脚本，则有一个此脚本中使用的“.to_crs（）”的已知错误。如果出现错误发生，您必须使用命令重置操作系统环境 ''os.environ[“PROJ_LIB”] = r“path-to-share-folder-in-environment”''.

合并多个文件

stack函数具有可选的输出参数，您可以在其中编写栅格添加到文件夹中的 TIFF 文件。如果要使用此功能，请确保有是要将 TIFF 文件写入的文件夹。 Stack 函数还返回两个对象，一个数组和一个 RasterIO 配置文件。做肯定会在变量中同时捕获。

In [14]:

# Stack Landsat bands
import glob
paths = glob.glob('/home/mw/input/precip7227/RainStormChina/RainStormChina/*/StormPeak*.tif')

out_path = './stacked_rasters.tif'  # 假设我们要输出的文件名为stacked_rasters.tif
array, raster_prof = es.stack(paths, out_path=out_path)

文件列表

In [15]:

paths

['/home/mw/input/precip7227/RainStormChina/RainStormChina/2014/StormPeak2014.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2005/StormPeak2005.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2017/StormPeak2017.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2007/StormPeak2007.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2015/StormPeak2015.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2006/StormPeak2006.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2013/StormPeak2013.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2002/StormPeak2002.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2012/StormPeak2012.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2003/StormPeak2003.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2011/StormPeak2011.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2004/StormPeak2004.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2001/StormPeak2001.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2018/StormPeak2018.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2019/StormPeak2019.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2010/StormPeak2010.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2008/StormPeak2008.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2009/StormPeak2009.tif',
 '/home/mw/input/precip7227/RainStormChina/RainStormChina/2016/StormPeak2016.tif']

创建范围对象

要获取栅格范围，plotting_extent请使用来自 ''es.stack（）'' 和 Rasterio 配置文件或元数据对象的数组。函数需要一个 Numpy 数组的层，这就是我们使用 ''arr[0]'' 的原因。该功能还需要 Rasterio 对象的空间变换，可以通过访问 Rasterio 配置文件中的 ''“transform”'' 键。

In [16]:

extent = plotting_extent(array[0], raster_prof["transform"])

绘制未裁剪的数据

您可以使用 ''ep.plot_rgb（）'' 查看裁剪前的边界和栅格请注意，数据似乎被洗掉了。

In [17]:

fig, ax = plt.subplots(figsize=(12, 12))
ep.plot_rgb(
    array,
    ax=ax,
    stretch=True,
    extent=extent,
    str_clip=0.5,
    title="RGB Image of Un-cropped Raster",
)
plt.show()

探索数据中的值范围

您可以使用 EarthPy ''hist（）'' 探索数据中找到的值范围功能。您是否注意到任何可能影响拉伸的极端值的图像？

In [18]:

ep.hist(array)
plt.show()

如何掩码

''es.stack（）'' 可以处理栅格中的 ''nodata'' 值。要使用此内容，请执行此操作参数，指定 ''nodata=''。这将遮罩包含的每个像素指定的“nodata”值。输出将是一个 numpy 掩码数组。

In [19]:

import glob
paths = glob.glob('/home/mw/input/precip7227/RainStormChina/RainStormChina/*/StormCount*.tif')

array_nodata, raster_prof_nodata = es.stack(paths, nodata=0)

# View hist of data with nodata values removed
ep.hist(
    array_nodata)
plt.show()

# Recreate extent object for the No Data array

extent_nodata = plotting_extent(
    array_nodata[0], raster_prof_nodata["transform"]
)

绘制未裁剪的数据（掩码后）

删除 nodata 值后，再次绘制数据。

In [20]:

fig, ax = plt.subplots(figsize=(12, 12))
ep.plot_rgb(
    array_nodata,
    ax=ax,
    stretch=True,
    extent=extent,
    str_clip=0.5,
    title="RGB image of Un-cropped Raster, No Data Value Selected",
)
plt.show()

裁剪数据

有时，您拥有的数据大于研究区域。在处理之前，先将数据裁剪到研究区域会更有效它在 Python 中。最快、最有效的选择是裁剪每个文件单独地将裁剪后的栅格写入新文件，然后堆叠将新文件放在一起。为此，请确保您具有 ShapeFile 边界以 GeoPandas 对象的形式，您可以用作裁剪对象。然后，循环浏览您要裁剪的每个文件并裁剪图像，然后将其写出到文件中。获取创建的栅格并将它们堆叠起来，就像您在前面的示例中堆叠了条带。

临时制作一个天津市shp

In [28]:

import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
# 读取全国地图数据
china_map = gpd.read_file('/home/mw/input/china1656/china_map/china_map/China_Province_2022.shp')
tj = china_map[china_map['省'] == '天津市']
tj = tj.rename(columns={'省': 'Province','省级码':'code','省类型':'type'})
tj.to_file('./tj.shp')

In [29]:

crop_bound = gpd.read_file(
    "/home/mw/project/tj.shp"
)
crop_bound

	Province	code	type	ENG_NAME	VAR_NAME	FIRST_GID	FIRST_TYPE	year	geometry
0	???	120000	???	Tianjin	Ti?n J?n	120000	Municipality	2022	POLYGON ((117.56937 40.19153, 117.56744 40.189...

重新投影数据

注意

如果你使用的是 Windows，请确保在此处设置你的环境！

如果数据位于不同的坐标中，裁剪功能将无法正常工作参考系统（CRS）。要解决此问题，请务必重新投影裁剪图层以匹配栅格数据的 CRS。要重投影数据，请先从栅格剖面中获取栅格的 CRS 对象。然后使用它使用 geopandas ''.to_crs'' 方法重新投影。

In [32]:

with rio.open(paths[0]) as raster_crs:
    crop_raster_profile = raster_crs.profile
    crop_bound_utm13N = crop_bound.to_crs(crop_raster_profile["crs"])
crop_raster_profile["crs"]

CRS.from_epsg(4326)

批量裁剪

当您需要裁剪和堆叠一组图像时，最有效的方法是先裁剪每个图像，然后堆叠它。 ''es.crop_all（）'' 是一种快速裁剪图像中所有波段的有效方法。该函数会将裁剪的栅格写入目录并返回文件路径列表，然后可以与 ''es.stack（）''。

In [42]:

output_dir = './'
band_paths_list = es.crop_all(
    paths, output_dir, crop_bound_utm13N, overwrite=True
)
print(glob.glob('/home/mw/project/*'))

堆叠所有数据

裁剪数据后，您就可以创建新堆叠了。

In [46]:

band_paths_list = glob.glob('/home/mw/project/*crop.tif')
cropped_array, array_raster_profile = es.stack(band_paths_list, nodata=-9999)
crop_extent = plotting_extent(
    cropped_array[0], array_raster_profile["transform"]
)

# Plotting the cropped image
# sphinx_gallery_thumbnail_number = 5
fig, ax = plt.subplots(figsize=(12, 6))
crop_bound_utm13N.boundary.plot(ax=ax, color="red", zorder=10)
ep.plot_rgb(
    cropped_array,
    ax=ax,
    stretch=True,
    extent=crop_extent,
    title="Cropped Raster and Fire Boundary",
)
plt.show()

裁剪单波段

如果你只需要裁剪一个光栅图像，你可以使用 EarthPy 的 ''es.crop_image（）'' 函数。此函数获取 Rasterio 对象并将其裁剪为提供的空间范围。

In [45]:

# Open Landsat image as a Rasterio object in order to crop it
stack_band_paths = glob.glob('/home/mw/project/*crop.tif')
with rio.open(stack_band_paths[0]) as src:
    single_cropped_image, single_cropped_meta = es.crop_image(
        src, crop_bound_utm13N
    )

# Create the extent object
single_crop_extent = plotting_extent(
    single_cropped_image[0], single_cropped_meta["transform"]
)

# Plot the newly cropped image
fig, ax = plt.subplots(figsize=(12, 6))
crop_bound_utm13N.boundary.plot(ax=ax, color="red", zorder=10)
ep.plot_bands(
    single_cropped_image,
    ax=ax,
    extent=single_crop_extent,
    title="Single Cropped Raster and Fire Boundary",
)
plt.show()