【Python科研】如何使用Python计算年度和季节的平均降水栅格数据并进行批量裁剪

目录

1.环境准备

2.设置文件路径

3.读取矢量数据

4.定义年份和季节

5.创建输出文件夹

6.裁剪栅格数据的函数

7.计算和保存年度平均降水数据

8.计算和保存季节平均降水数据

9.结论

10.完整代码

        本次分享内容中，我们将介绍如何使用Python计算和裁剪年度和季节的降水栅格数据。本教程主要涉及到以下内容：

  1. 计算2000-2021年的年度和季节平均降水量；
  2. 使用研究区域矢量数据裁剪出研究区域的降水栅格数据。

        我们将使用 rasterio 和 geopandas 两个库来实现这些操作。

 图1|现有需要计算处理和裁剪的降水栅格数据

1.环境准备

        首先，确保你的Python环境中已经安装了必要的库。如果没有的话，你可以使用以下命令来安装这些库（在Anaconda进行）：

pip install rasterio geopandas numpy

2.设置文件路径

        将文件路径设置为你自己的数据文件夹路径和矢量数据路径：

import os
import rasterio
import numpy as np
import geopandas as gpd
from rasterio.mask import mask

input_folder = r"你的降水栅格数据文件夹路径"
output_folder = r"保存结果的文件夹路径"
clipped_output_folder = r"保存裁剪后结果的文件夹路径"
shapefile_path = r"矢量数据文件路径"

3.读取矢量数据

        使用 geopandas 读取研究区域的矢量数据：

shapefile = gpd.read_file(shapefile_path)

4.定义年份和季节

        设置要处理的年份和季节：

years = list(range(2000, 2021))  # 从2000到2020年
seasons = {
    "spring": [3, 4, 5],
    "summer": [6, 7, 8],
    "autumn": [9, 10, 11],
    "winter": [12, 1, 2]
}

5.创建输出文件夹

        确保输出文件夹存在：

os.makedirs(output_folder, exist_ok=True)
os.makedirs(clipped_output_folder, exist_ok=True)

6.裁剪栅格数据的函数

        定义一个函数来裁剪栅格数据，并设置无效值：

def clip_raster(raster, shapes, nodata_value):
    out_image, out_transform = mask(raster, shapes, crop=True, nodata=nodata_value, filled=False)
    out_meta = raster.meta.copy()
    out_meta.update({
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform,
        "nodata": nodata_value
    })
    return out_image, out_meta

7.计算和保存年度平均降水数据

逐年计算并保存年度平均降水数据，然后裁剪结果：

for year in years:
    annual_data = []
    for month in range(1, 13):
        file = f"{input_folder}/pre{year}{str(month).zfill(2)}.tif"
        if os.path.exists(file):
            with rasterio.open(file) as src:
                data = src.read(1).astype(np.float32)
                data[data == src.nodata] = np.nan  # 替换无效值为NaN
                annual_data.append(data)

    if annual_data:
        annual_mean = np.nanmean(np.array(annual_data), axis=0)  # 忽略NaN值
        
        # 保存年度平均降水数据
        output_file = f"{output_folder}/annual_mean_{year}.tif"
        with rasterio.open(
            output_file, 'w',
            driver='GTiff',
            height=annual_mean.shape[0],
            width=annual_mean.shape[1],
            count=1,
            dtype=annual_mean.dtype,
            crs=src.crs,
            transform=src.transform,
            nodata=np.nan
        ) as dst:
            dst.write(annual_mean, 1)
        
        # 裁剪年度平均降水数据
        with rasterio.open(output_file) as src:
            out_image, out_meta = clip_raster(src, shapefile.geometry, nodata_value=np.nan)
            clipped_output_file = f"{clipped_output_folder}/annual_mean_{year}_clipped.tif"
            with rasterio.open(clipped_output_file, 'w', **out_meta) as dst:
                dst.write(out_image[0], 1)

8.计算和保存季节平均降水数据

        逐年逐季计算并保存季节平均降水数据，然后裁剪结果：

for year in years:
    for season, months in seasons.items():
        season_data = []
        for month in months:
            if month == 1 or month == 2:  # 处理跨年的冬季月份
                season_year = year + 1
            else:
                season_year = year
            
            file = f"{input_folder}/pre{season_year}{str(month).zfill(2)}.tif"
            if os.path.exists(file):
                with rasterio.open(file) as src:
                    data = src.read(1).astype(np.float32)
                    data[data == src.nodata] = np.nan  # 替换无效值为NaN
                    season_data.append(data)

        if season_data:
            season_mean = np.nanmean(np.array(season_data), axis=0)  # 忽略NaN值
            
            # 保存季节平均降水数据
            output_file = f"{output_folder}/season_mean_{season}_{year}.tif"
            with rasterio.open(
                output_file, 'w',
                driver='GTiff',
                height=season_mean.shape[0],
                width=season_mean.shape[1],
                count=1,
                dtype=season_mean.dtype,
                crs=src.crs,
                transform=src.transform,
                nodata=np.nan
            ) as dst:
                dst.write(season_mean, 1)
            
            # 裁剪季节平均降水数据
            with rasterio.open(output_file) as src:
                out_image, out_meta = clip_raster(src, shapefile.geometry, nodata_value=np.nan)
                clipped_output_file = f"{clipped_output_folder}/season_mean_{season}_{year}_clipped.tif"
                with rasterio.open(clipped_output_file, 'w', **out_meta) as dst:
                    dst.write(out_image[0], 1)

9.结论

        通过上述步骤，我们可以计算年度和季节的平均降水量，并使用矢量数据裁剪结果，使其仅包含研究区域内的有效像素，这样可以更准确地分析研究区域的降水情况。

10.完整代码

import os
import rasterio
import numpy as np
import geopandas as gpd
from rasterio.mask import mask

# 设置文件路径
input_folder = r"你的降水栅格数据文件夹路径"
output_folder = r"保存结果的文件夹路径"
clipped_output_folder = r"保存裁剪后结果的文件夹路径"
shapefile_path = r"矢量数据文件路径"

# 读取矢量数据
shapefile = gpd.read_file(shapefile_path)

# 定义年份和季节
years = list(range(2000, 2021))  # 从2000到2020年
seasons = {
    "spring": [3, 4, 5],
    "summer": [6, 7, 8],
    "autumn": [9, 10, 11],
    "winter": [12, 1, 2]
}

# 创建输出文件夹（如果不存在）
os.makedirs(output_folder, exist_ok=True)
os.makedirs(clipped_output_folder, exist_ok=True)

# 裁剪栅格数据的函数
def clip_raster(raster, shapes, nodata_value):
    out_image, out_transform = mask(raster, shapes, crop=True, nodata=nodata_value, filled=False)
    out_meta = raster.meta.copy()
    out_meta.update({
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform,
        "nodata": nodata_value
    })
    return out_image, out_meta

# 计算年度平均降水数据
for year in years:
    annual_data = []
    for month in range(1, 13):
        file = f"{input_folder}/pre{year}{str(month).zfill(2)}.tif"
        if os.path.exists(file):
            with rasterio.open(file) as src:
                data = src.read(1).astype(np.float32)
                data[data == src.nodata] = np.nan  # 替换无效值为NaN
                annual_data.append(data)

    if annual_data:
        annual_mean = np.nanmean(np.array(annual_data), axis=0)  # 忽略NaN值
        
        # 保存年度平均降水数据
        output_file = f"{output_folder}/annual_mean_{year}.tif"
        with rasterio.open(
            output_file, 'w',
            driver='GTiff',
            height=annual_mean.shape[0],
            width=annual_mean.shape[1],
            count=1,
            dtype=annual_mean.dtype,
            crs=src.crs,
            transform=src.transform,
            nodata=np.nan
        ) as dst:
            dst.write(annual_mean, 1)
        
        # 裁剪年度平均降水数据
        with rasterio.open(output_file) as src:
            out_image, out_meta = clip_raster(src, shapefile.geometry, nodata_value=np.nan)
            clipped_output_file = f"{clipped_output_folder}/annual_mean_{year}_clipped.tif"
            with rasterio.open(clipped_output_file, 'w', **out_meta) as dst:
                dst.write(out_image[0], 1)

# 计算季节平均降水数据
for year in years:
    for season, months in seasons.items():
        season_data = []
        for month in months:
            if month == 1 or month == 2:  # 处理跨年的冬季月份
                season_year = year + 1
            else:
                season_year = year
            
            file = f"{input_folder}/pre{season_year}{str(month).zfill(2)}.tif"
            if os.path.exists(file):
                with rasterio.open(file) as src:
                    data = src.read(1).astype(np.float32)
                    data[data == src.nodata] = np.nan  # 替换无效值为NaN
                    season_data.append(data)

        if season_data:
            season_mean = np.nanmean(np.array(season_data), axis=0)  # 忽略NaN值
            
            # 保存季节平均降水数据
            output_file = f"{output_folder}/season_mean_{season}_{year}.tif"
            with rasterio.open(
                output_file, 'w',
                driver='GTiff',
                height=season_mean.shape[0],
                width=season_mean.shape[1],
                count=1,
                dtype=season_mean.dtype,
                crs=src.crs,
                transform=src.transform,
                nodata=np.nan
            ) as dst:
                dst.write(season_mean, 1)
            
            # 裁剪季节平均降水数据
            with rasterio.open(output_file) as src:
                out_image, out_meta = clip_raster(src, shapefile.geometry, nodata_value=np.nan)
                clipped_output_file = f"{clipped_output_folder}/season_mean_{season}_{year}_clipped.tif"
                with rasterio.open(clipped_output_file, 'w', **out_meta) as dst:
                    dst.write(out_image[0], 1)

print("年度和季节平均降水数据计算和裁剪完成！")

图2|计算裁剪之后得到的研究区域降水数据

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：/a/726664.html

如若内容造成侵权/违法违规/事实不符，请联系我们进行投诉反馈qq邮箱809451989@qq.com，一经查实，立即删除！

相关文章

复合构件之消息对话框

复合构件之消息对话框

代码; #include <gtk-2.0/gtk/gtk.h> static void font_dialog_response(GtkFontSelectionDialog *dialog,gint response,gpointer data)// 处理字体选择对话框按钮按下事件 {gchar *font;GtkWidget *message;switch (response) {case (GTK_RESPONSE_APPLY):case (GTK_R…
阅读更多...
ClipArt ETC - 典雅的剪贴画网站

ClipArt ETC - 典雅的剪贴画网站

文章目录 ClipArt ETCClippix佛罗里达教学技术中心课堂数字内容 ClipArt ETC 网站地址&#xff1a; https://etc.usf.edu/clipart/ ClipArt ETC为学生和教师提供了超过71,500件高质量的教育剪贴画。 每个插图都有图像大小的选择以及学校项目中正确引用的完整源信息。 所有图像…
阅读更多...
UniApp 开发微信小程序教程（一）：准备工作和环境搭建，项目结构和配置

UniApp 开发微信小程序教程（一）：准备工作和环境搭建，项目结构和配置

文章目录 一、准备工作和环境搭建1. 安装 HBuilderX步骤&#xff1a; 2. 注册微信开发者账号步骤&#xff1a; 3. 创建 UniApp 项目步骤&#xff1a; 二、项目结构和配置1. UniApp 项目结构2. 配置微信小程序修改 manifest.json修改 pages.json 3. 添加首页文件index.vue 示例&…
阅读更多...
【查缺补漏】python

【查缺补漏】python

python查缺补漏 底板除 还有一种除法是//&#xff0c;称为地板除&#xff0c;两个整数的除法仍然是整数&#xff1a; >>> 10 // 3 3你没有看错&#xff0c;整数的地板除//永远是整数&#xff0c;即使除不尽。要做精确的除法&#xff0c;使用/就可以。 因为//除法只…
阅读更多...
一文讲清楚分销裂变是什么？怎么做好分销裂变？【附案例】

一文讲清楚分销裂变是什么？怎么做好分销裂变？【附案例】

在数字化营销日益盛行的今天&#xff0c;分销裂变作为一种高效的推广手段&#xff0c;受到了越来越多企业的青睐。那么&#xff0c;分销裂变究竟是什么&#xff1f;我们又该如何做好分销裂变呢&#xff1f;林叔将从定义、方法以及案例分析三个方面进行阐述。 一、分销裂变是什…
阅读更多...
MySQL的数据存储一定是基于硬盘吗？

MySQL的数据存储一定是基于硬盘吗？

一、典型回答 不是的&#xff0c;MySQL也可以基于内存的&#xff0c;即MySQL的内存表技术。它允许将数据和索引存储在内存中&#xff0c;从而提高了检验速度和修改数据的效率。优点包括具有快速响应的查询性能和节约硬盘存储空间。此外&#xff0c;使用内存表还可以实现更高的复…
阅读更多...
数据库讲解---(数据库设计)

数据库讲解---(数据库设计)

目录 一.数据库设计概述 1.1数据库设计的内容 1.1.1数据库的结构设计 1.1.2数据库的行为设计 1.2数据库设计方法 1.2.1直观设计法 1.2.2规范设计法 1.2.3计算机辅助设计法 1.2.4自动化设计法 1.3数据库设计的基本步骤 1.3.1需求分析 1.3.2概念结构设计 1.3.3逻辑结…
阅读更多...
丹尼尔·T·琼斯：精益生产到底是什么？

丹尼尔·T·琼斯：精益生产到底是什么？

本文摘要自《精益思想》、《改变世界的机器》作者之一丹尼尔T琼斯的文章。丹尼尔T琼斯是一位学者、英国作家和研究员。他曾多次获得瑞士山吉奥卓越运营奖研究与专业出版类别的奖项&#xff0c;也包括了国际精益六西格玛研究所&#xff08;ILSSI&#xff09;[1]的"精益思想…
阅读更多...
CentOS Linux 7系统中离线安装MySQL5.7步骤

CentOS Linux 7系统中离线安装MySQL5.7步骤

预计数据文件存储目录为&#xff1a;/opt/mysql/data 1、文件下载&#xff1a; 安装文件下载链接&#xff1a;https://downloads.mysql.com/archives/community/ 2、检查当前系统是否安装过MySQL [rootcnic51 mysql]# rpm -qa|grep mariadb mariadb-libs-5.5.68-1.el7.x86_6…
阅读更多...
Java中的运算符及其示例

Java中的运算符及其示例

Java中的运算符及其示例 运算符是指示编译器执行特定操作的符号。例如&#xff0c;“”运算符指示编译器执行加法&#xff0c;“>”运算符指示编译执行比较&#xff0c;“”用于赋值等等。在本指南中&#xff0c;我们将借助示例讨论java中的操作。 运算符和操作数&#…
阅读更多...
一文读懂什么是SaaS产品运营？如何做好SaaS产品运营？

一文读懂什么是SaaS产品运营？如何做好SaaS产品运营？

在当今数字化时代&#xff0c;SaaS&#xff08;Software-as-a-Service&#xff0c;软件即服务&#xff09;产品已成为企业运营不可或缺的一部分。本文将结合具体案例&#xff0c;深入解析SaaS产品运营的定义与策略。 一、什么是SaaS产品运营&#xff1f; SaaS产品运营是指通过…
阅读更多...
由于“xinput1_3.dll缺失“而导致的错误有哪些解决办？分享几种修复xinput1_3.dll丢失的方法

由于“xinput1_3.dll缺失“而导致的错误有哪些解决办？分享几种修复xinput1_3.dll丢失的方法

当您尝试运行某些游戏或程序时&#xff0c;可能会遇到由于"xinput1_3.dll缺失"而导致的错误。这个DLL文件是MicrosoftDirectX的一部分&#xff0c;用于处理游戏中的输入设备&#xff0c;如操纵杆和游戏手柄。下面我们将探讨为何电脑会缺少xinput1_3.dll文件&#xff…
阅读更多...
Spring系统学习 - FactoryBean和基于XML的自动装配

Spring系统学习 - FactoryBean和基于XML的自动装配

Factory Bean Spring的FactoryBean是一个特殊的Bean&#xff0c;用于创建其他Bean实例。FactoryBean接口定义了一个工厂Bean&#xff0c;该Bean可以用来生成其他Bean的实例。通过实现FactoryBean接口&#xff0c;开发人员可以自定义Bean的创建逻辑&#xff0c;实现更灵活的Bea…
阅读更多...
软硬件节水“组合拳”,助力智慧灌区信息化建设！

软硬件节水“组合拳”,助力智慧灌区信息化建设！

水资源短缺已成为全球共同面临的挑战&#xff0c;尤其在农业灌溉领域&#xff0c;其影响尤为显著。农业作为水资源消耗的主要行业之一&#xff0c;在日益严峻的水资源形势下&#xff0c;构建节水型灌区的紧迫性日益凸显。 节水型灌区的建设&#xff0c;旨在通过优化灌溉方式、减…
阅读更多...
【C++ | 移动构造函数】C++11的 移动构造函数 详解及例子代码

【C++ | 移动构造函数】C++11的 移动构造函数 详解及例子代码

&#x1f601;博客主页&#x1f601;&#xff1a;&#x1f680;https://blog.csdn.net/wkd_007&#x1f680; &#x1f911;博客内容&#x1f911;&#xff1a;&#x1f36d;嵌入式开发、Linux、C语言、C、数据结构、音视频&#x1f36d; ⏰发布时间⏰&#xff1a;2024-06-12 2…
阅读更多...
「Python-docx 专栏」docx 获取页面大小、设置页面大小（纸张大小）

「Python-docx 专栏」docx 获取页面大小、设置页面大小（纸张大小）

本文目录 前言一、docx纸张大小介绍1、document.xml① 关于 document.xml 的一些知识点② 纸张大小在哪里③ 纸张大小都有啥④ EMU对应的尺寸列表二、获取docx纸张大小1、完整代码2、运行效果图三、python为docx设置纸张大小1、完整代码2、效果图前言 今天的这边文章，我们来说…
阅读更多...
Java多线程基础知识-2

Java多线程基础知识-2

线程的3个方法&#xff1a; Thread.sleep()&#xff1a;当前线程睡眠多少毫秒&#xff0c;让给其他线程去执行。 Thread.yield()&#xff1a;当前线程退出一下&#xff0c;进入到等待队列&#xff0c;让其他线程执行&#xff0c;即让出线程一下。 Thread.join()&#xff1a;…
阅读更多...
LabVIEW_TDMS

LabVIEW_TDMS

1.TDMS设置属性 想给这里写属性怎么整 使用TDMS设置属性函数时&#xff0c;对组名称与通道名称不设置&#xff0c;即可达到上图中的样式。 PS&#xff1a;属性名称如果设置一样则最终生效的值为最后写入的值。如将属性2修改为属性1&#xff0c;则最终只有1个属性1&#xff0c…
阅读更多...
Pikachu靶场--文件上传

Pikachu靶场--文件上传

参考借鉴 Pikachu靶场之文件上传漏洞详解_皮卡丘文件上传漏洞-CSDN博客 文件上传漏洞&#xff1a;pikachu靶场中的文件上传漏洞通关_pikachu文件上传通关-CSDN博客 client check 在桌面新建一个文件夹&#xff0c;准备一个hello.php文件&#xff0c;文件写入如下代码 <?p…
阅读更多...
Safari浏览器下载文件时，文件名会URL encoded

Safari浏览器下载文件时，文件名会URL encoded

问题&#xff1a;相同链接下载文件&#xff0c;safari文件名编码异常 解决&#xff1a;response.setHeader("Content-Disposition", "attachment;filename*utf-8" URLEncoder.encode(filename, "UTF-8")); 问题描述 谷歌下载&#xff08;正常&a…
阅读更多...
最新文章

Copyright @ 2022~2023