Loading [MathJax]/jax/output/CommonHTML/config.js
首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >专栏 >使用 EarthPy 堆叠和裁剪tif栅格数据

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

作者头像
用户11172986
发布于 2024-06-20 11:17:34
发布于 2024-06-20 11:17:34
24600
代码可运行
举报
文章被收录于专栏:气python风雨气python风雨
运行总次数:0
代码可运行

使用 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 [ ]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
!pip install earthpy -i https://pypi.mirrors.ustc.edu.cn/simple/
代码语言:javascript
代码运行次数:0
运行
复制

In [1]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
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
代码语言:javascript
代码运行次数:0
运行
复制

注意

注意

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

合并多个文件

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

In [14]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
# 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)
代码语言:javascript
代码运行次数:0
运行
复制

文件列表

In [15]:

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
paths
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
['/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]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
extent = plotting_extent(array[0], raster_prof["transform"])
代码语言:javascript
代码运行次数:0
运行
复制

绘制未裁剪的数据

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

In [17]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
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()
代码语言:javascript
代码运行次数:0
运行
复制

探索数据中的值范围

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

In [18]:

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
ep.hist(array)
plt.show()

如何掩码

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

In [19]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
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"]
)
代码语言:javascript
代码运行次数:0
运行
复制

绘制未裁剪的数据(掩码后)

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

In [20]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
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()
代码语言:javascript
代码运行次数:0
运行
复制

裁剪数据

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

临时制作一个天津市shp

In [28]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
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')
代码语言:javascript
代码运行次数:0
运行
复制

In [29]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
crop_bound = gpd.read_file(
    "/home/mw/project/tj.shp"
)
crop_bound
代码语言:javascript
代码运行次数:0
运行
复制

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]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
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"]
代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
CRS.from_epsg(4326)

批量裁剪

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

In [42]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
output_dir = './'
band_paths_list = es.crop_all(
    paths, output_dir, crop_bound_utm13N, overwrite=True
)
print(glob.glob('/home/mw/project/*'))
代码语言:javascript
代码运行次数:0
运行
复制

堆叠所有数据

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

In [46]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
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()
代码语言:javascript
代码运行次数:0
运行
复制

裁剪单波段

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

In [45]:

代码语言:javascript
代码运行次数:0
运行
复制
代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
# 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()
代码语言:javascript
代码运行次数:0
运行
复制

小结

以上是基于earthpy的官方教程修改,因为官方的示例数据较难下载就改成自己的数据

看得出erathpy对于tif数据处理较为便利,堆叠和掩膜都比较简短

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2024-05-16,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 气python风雨 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
暂无评论
推荐阅读
编辑精选文章
换一批
多年暴雨tif数据集合成为一个nc数据
当处理多年暴雨的 TIF 数据集时,我们可以使用 rioxarray 库将这些数据合成为一个 NetCDF (nc) 文件。NetCDF 是一种常用的科学数据格式,它具有跨平台、可扩展和自描述的特点,非常适合存储和共享地理空间数据。
用户11172986
2024/06/20
4490
多年暴雨tif数据集合成为一个nc数据
进阶!dask解决超高精度tif读取与绘图难问题
又有读者来信 要求如下: 希望小编帮忙看看能不能解决。是关于能不能在已经截取出来的省份中添加对应的dem地形呢,并且根据需要添加上需要的城市所在的地理位置,比如在已绘制的图中标注出三亚的所在地
用户11172986
2024/06/20
3110
进阶!dask解决超高精度tif读取与绘图难问题
两个气象雷达库的简单使用
项目地址:https://pycwr.readthedocs.io/en/latest/draw.html
用户11172986
2024/06/20
6270
两个气象雷达库的简单使用
Earthpy | 这样超赞的艺术地图也能轻松绘制...
最近在整理Python数据可视化课程的拓展内容时,发现了一个处理空间数据的超赞工具-「earthpy」,也解决了一个绘制艺术地图的问题,下面就给大家详细介绍一下这个工具~~
DataCharm
2024/05/11
3420
Earthpy | 这样超赞的艺术地图也能轻松绘制...
台风路径可视化 | 强台风“天兔”移动轨迹图(Chiikawa版)
大家好!最近,2024年第25号台风天兔引起了网友们的广泛关注,因为它有一个特别的名字——“ウサギ”(Usagi),直译过来就是“兔子”。其本义是指天兔星座。 由于这个名字与《Chiikawa》漫画中的可爱角色“乌萨奇”(Usagi)同名,网友们亲切地称其为“乌萨奇台风”。这不仅增加了台风话题的趣味性,也让许多人对这次台风产生了更多的关注。
用户11172986
2024/11/25
6830
台风路径可视化 | 强台风“天兔”移动轨迹图(Chiikawa版)
GPM 卫星三级产品的简单可视化
之前GPM数据写了下载教程,现在简单试试可视化,毕竟是nc格式数据(下载可选),用起来相对简单
用户11172986
2024/06/20
2260
GPM 卫星三级产品的简单可视化
简单实现区域shp合并
偶然看见有人在求什么西南区域,东北区域的shp,写一期不求人攻略。 前面写过怎么裁剪,这次讲讲怎么合并,实现区域shp自由
用户11172986
2024/06/20
2130
简单实现区域shp合并
葵花八号AHI真彩图可视化
📢 版权声明:公益性质转载需联系作者本人获取授权。转载本文时,请务必文字注明“来自:和鲸社区:酷炫用户名”,并附带本项目超链接。
用户11172986
2024/07/18
4041
葵花八号AHI真彩图可视化
metpy绘制锋生与冷锋
https://www.heywhale.com/mw/project/65485a22d74b63fed5f03f49
用户11172986
2024/06/20
4640
metpy绘制锋生与冷锋
metpy函数平滑台风风场流线图
九点平滑的工作原理是将风速数据中的每个值替换为该值及其八个相邻值的平均值。这具有平滑数据和消除任何高频噪声的效果。 下面是一步一步解释九点平滑器是如何工作的: 创建一个新数组来存储平滑后的值。 迭代风速数据数组。 对于数组中的每个值,获取八个相邻值。 计算九个值的平均值。 在新数组中存储平均值。 平滑过程完成后,新阵列将包含平滑后的风速数据。
用户11172986
2024/06/20
2180
metpy函数平滑台风风场流线图
用R处理NASA数据(.hdf 或.nc文件)
这里不在赘述,参考如何获取NASA数据,下面的例子根据下载的LandCover与Rainfall数据进行展示,如何利用R语音进行读取,然后绘图。先加载所需R包及地图文件
Jamesjin63
2022/10/25
1.4K0
用R处理NASA数据(.hdf 或.nc文件)
GPM卫星数据hdf5格式读取与绘图
你刚开始拿到数据多半不知怎么看结构,一定很疑惑f['Swath/latentHeating'][:]怎么来的 hdf5数据逻辑和nc不太一样, 且看我下面如何操作
用户11172986
2024/06/20
6180
GPM卫星数据hdf5格式读取与绘图
地科Python数据分析案例 | 绘制黄土高原局部区域的沟壑覆盖度分析图
黄土高原位于中国中部偏北部,为中国四大高原之一。黄土高原是世界上水土流失最严重和生态环境最脆弱的地区之一,除许多石质山地外,大部分区域为厚层黄土覆盖,经流水长期强烈侵蚀,逐渐形成千沟万壑、地形支离破碎的特殊自然景观。
陈南GISer
2023/11/03
1.2K0
地科Python数据分析案例 | 绘制黄土高原局部区域的沟壑覆盖度分析图
两种micaps站点数据的简单绘制方法
用户11172986
2024/06/20
4600
两种micaps站点数据的简单绘制方法
Python绘制地图专用库(Cartopy)
地图绘制 大家在绘制栅格地图的时候有可能还在使用ArcGIS进行出图,但是ArcGIS出图比较慢,而且批量出图的时候又比较麻烦。 今天给大家介绍一个Python中用于地图绘制的库,Cartopy,这个库跟basemap非常相似,不过basemap现在已经不再更新。所以大家使用Python绘制地图还是使用Cartopy比较好。 Cartopy简介 Cartopy是一个Python软件包,用于地理空间数据处理,以便生成地图和其他地理空间数据分析。 网址:https://scitools.org.uk/carto
GIS与遥感开发平台
2022/04/29
2.8K0
Python绘制地图专用库(Cartopy)
meteva站点插值填色与白化
说到插值大概会提到日常用的scipy的linear和cubic,克里金插值等等 meteva也有插值功能,不论是站点插网格,网格插站点,还是网格插网格统统都有 本文主要测试meteva的IDW与cressman站点插值 并基于插值后的数据测试插值后的白化效果
用户11172986
2024/06/20
2940
meteva站点插值填色与白化
常见地图白化方法(一)
地图白化是一种绘制地图的技术,它可以实现对感兴趣区域以外的数据进行遮盖或填充白色的效果,从而突出显示目标区域的特征。 地图白化的原理是利用 shapefile 文件中的多边形坐标来创建一个剪切路径,然后将这个路径应用到 matplotlib 的绘图对象上,使得只有路径内的数据可见,路径外的数据被隐藏或覆盖。 气象家园的五星上将兰溪之水说过,其实所谓的“白化”一般就两种途径:①mask数据;②mask图形
用户11172986
2024/06/20
2540
常见地图白化方法(一)
一瞬又一瞬,累积起来便是一生 | ERA5数据计算垂直积分整层水汽通量散度
由于可视化代码过长隐藏,可点击运行Fork查看 若没有成功加载可视化图,点击运行可以查看 ps:隐藏代码在【代码已被隐藏】所在行,点击所在行,可以看到该行的最右角,会出现个三角形,点击查看即可
用户11172986
2024/12/30
8780
一瞬又一瞬,累积起来便是一生 | ERA5数据计算垂直积分整层水汽通量散度
WRFOUT 涡度平流和温度平流计算与可视化
涡度平流和温度平流是两种常见的气象诊断量,可以帮助我们更好地理解大气运动和热力学过程。 以下代码将计算上述气象诊断量并可视化。
用户11172986
2024/06/20
6940
WRFOUT 涡度平流和温度平流计算与可视化
求栅格序列每个像元的变化趋势和对应P值
讲完了geotiff格式数据的读取和保存,本文讲下怎么用python处理一系列的栅格数据(本文以时间序列为例)。
自学气象人
2022/11/14
3.1K0
求栅格序列每个像元的变化趋势和对应P值
相关推荐
多年暴雨tif数据集合成为一个nc数据
更多 >
交个朋友
加入腾讯云官网粉丝站
蹲全网底价单品 享第一手活动信息
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档
本文部分代码块支持一键运行,欢迎体验
本文部分代码块支持一键运行,欢迎体验