python--读取TRMM-3B43月平均降水绘制气候态空间分布图（陆地区域做掩膜）

news2025/1/13 16:52:55

python–读取TRMM-3B43月平均降水绘制气候态空间分布图（陆地区域做掩膜）

成果展示

在这里插入图片描述

TRMM降水数据介绍

热带降雨测量任务(The Tropical Rainfall Measuring Mission，TRMM)是美国国家航空航天局(NASA)和日本国家太空发展署(National Space Development Agency)的一项联合太空任务，旨在监测和研究热带和亚热带降雨及其相关的能量释放。该任务使用5种仪器: 降水雷达(Precipitation Radar，PR)、 TRMM 微波成像仪(TRMM Microwave Imager，TMI)、可见红外扫描仪(Visible Infrared Scanner，VIRS)、云与地球辐射能量系统(Clouds & Earths Radiant Energy System，CERES)和闪电成像传感器(Lightning Imaging Sensor，LSI)。TMI 和 PR 是用于降水的主要仪器。这些仪器被用于形成 TRMM 多卫星降水分析(TRMM Multi-satellite Precipitation Analysis，TMPA)的 TRMM 组合仪器(TRMM Combined Instrument ，TCI)校准数据集(TRMM 2B31)的算法中，其 TMPA 3B43月平均降水量和 TMPA 3B42日平均和次日(3小时)平均是最相关的 TRMM 相关气候研究的产品。3B42和3B43的空间分辨率为0.25 ° ，1998年至今覆盖北纬50 ° 至南纬50 ° 。

本文中用到的数据主要为TRMM-3B43月平均产品，用于绘制降水的气候态空间分布图

TRMM-3B43产品如下所示：

  • 空间分辨率：0.25°
  • 时间覆盖范围：1999.01 - 2020.01
  • 经纬度范围： 经度：0-360°，纬度：-50°S-50°N
  • 单位为：mm/hour
  • 降水类型 ： 累计降水

这里使用的python系统环境：linux系统，因为windows上好多库都不好用

关于掩膜的方法，之前出过教程了，这里不再重复：

python 对陆地数据进行掩膜的两种方法

在这里插入图片描述
在这里插入图片描述

官网数据下载链接

数据处理

  • 这里所使用的3B43降水资料数据为.HDF格式，因此需要使用pyhdf这个库来读取，不习惯的可以下载netcdf的数据格式
  • 由于数据中缺少经纬度信息（也可能是我没有找到），为了实现区域切片，这里手动造了一个dataarray的数据，从而实现切片的过程
  • 数据中读取的变量为precipitation，读取完之后是个二维的数组，为了给他加上时间纬度，所以手动给他进行了扩维，之后实现多年的气候态平均计算
  • 使用global_land_mask实现对于陆地区域的掩膜处理
  • 使用cnmaps 实现中国的区域绘制，这里的cnmaps的库在windows上可能不好安装，也直接使用cartopyax.coastlines('50m')自带的海岸线
  • 将原始数据单位转化为 mm/year，这里只是简单的转换*24*365

代码实现

1、首先读取10年的数据路径

import xarray as xr
import os,glob
import numpy as np
from pyhdf.SD import SD
path = '/Datadisk/TRMM/3B43/'
file_list = []
for year in range(2009,2019):
    folder = os.path.join(path, str(year))
    file_name = glob.glob(folder+'/3B43.'+str(year)+'*.HDF')
    file_name.sort()
    file_list.extend(file_name) 
file_list.sort()

在这里插入图片描述
2、封装数据读取函数，并对需要的区域进行切片，我选择的区域为经度：[90.0, 145]，纬度：[-10, 55]，并将循环读取的10年月平均数组创建为dataarray的格式方面后续掩膜，这里的时间可以通过pandas自己创建时间序列，我这里偷懒直接读取了之前处理过的月平均gpcp的time了

dt = xr.open_dataset("/gpcp_monthly_mask.nc")
precip = dt.sel(time=slice('2009','2018'),lat =slice(-10,55),lon = slice(90,145)).precip
time_num = precip.time.values
def get_data(path,time):
    da = SD(path)
    pre = da.select('precipitation')[:]
    pre = np.expand_dims(pre,axis=0)
    lon = np.arange(-180,180.,0.25)
    lat = np.arange(-50,50.,0.25)
    time = time
#   time = datetime.datetime.utcfromtimestamp(tim).strftime('%Y-%m-%d %H:%M:%S')
    da = xr.DataArray(pre, 
                          dims=['time','lon','lat'],
                          coords=dict(
                              lon=(['lon'], lon),
                              lat=(['lat'], lat),
                              time=(['time'],[time])),
                              )
    ##############################################################################
    lon_range = [90.0, 145]
    lat_range = [-10, 55]
    da = da.sel(lon=slice(*lon_range), lat=slice(*lat_range))
    x ,y = da.lon,da.lat   
    return da,x,y
rain = np.zeros((len(file_list),221,240))            
for i in range(len(file_list)):
    print(i)
    da,x,y =  get_data(file_list[i], time_num[i])
    rain[i] = da
ds = xr.DataArray(rain, 
                      dims=['time','lon','lat'],
                      coords=dict(
                          lon=(['lon'], x.data),
                          lat=(['lat'], y.data),
                          time=(['time'],time_num)),
                          )

3、对数据的陆地部分进行掩膜，并计算气候态平均，转换单位为mm/year

from global_land_mask import globe
def mask_land(ds, label='land', lonname='lon'):
    if lonname == 'lon':
        lat = ds.lat.data
        lon = ds.lon.data
        if np.any(lon > 180):
            lon = lon - 180
            lons, lats = np.meshgrid(lon, lat)
            mask = globe.is_ocean(lats, lons)
            temp = []
            temp = mask[:, 0:(len(lon) // 2)].copy()
            mask[:, 0:(len(lon) // 2)] = mask[:, (len(lon) // 2):]
            mask[:, (len(lon) // 2):] = temp
        else:
            lons, lats = np.meshgrid(lon, lat)# Make a grid
            mask = globe.is_ocean(lats, lons)# Get whether the points are on ocean.
        ds.coords['mask'] = (('lat', 'lon'), mask)
    elif lonname == 'longitude':
        lat = ds.latitude.data
        lon = ds.longitude.data
        if np.any(lon > 180):
            lon = lon - 180
            lons, lats = np.meshgrid(lon, lat)
            mask = globe.is_ocean(lats, lons)
            temp = []
            temp = mask[:, 0:(len(lon) // 2)].copy()
            mask[:, 0:(len(lon) // 2)] = mask[:, (len(lon) // 2):]
            mask[:, (len(lon) // 2):] = temp
        else:
            lons, lats = np.meshgrid(lon, lat)
            mask = globe.is_ocean(lats, lons)
        lons, lats = np.meshgrid(lon, lat)
        mask = globe.is_ocean(lats, lons)
        ds.coords['mask'] = (('latitude', 'longitude'), mask)
    if label == 'land':
        ds = ds.where(ds.mask == True)
    elif label == 'ocean':
        ds = ds.where(ds.mask == False)
    return ds
data =  mask_land(ds,'land')
precip_mean = np.nanmean(data*24*365,axis=0)

4、绘图，保存图片

import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import cmaps
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cnmaps import get_adm_maps, draw_maps
box = [100,140,0,50]
xstep,ystep = 10,10
proj = ccrs.PlateCarree(central_longitude=180)
plt.rcParams['font.family'] = 'Times New Roman',
fig = plt.figure(figsize=(8,7),dpi=200)
fig.tight_layout()
ax = fig.add_axes([0.1,0.2,0.8,0.7],projection = proj)
ax.set_extent(box,crs=ccrs.PlateCarree())
draw_maps(get_adm_maps(level='国')) #这里如果库不好安装的话可以使用下面注释的代码，cartopy自带的海岸线
# ax.coastlines('50m')
ax.set_xticks(np.arange(box[0],box[1]+xstep, xstep),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(box[2], box[3]+1, ystep),crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=False)#True/False
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_title('TRMM(mm/year)',fontsize=16,pad=8,loc='left')
ax.tick_params(    which='both',direction='in',
                width=0.7,
                    pad=8, 
                    labelsize=14,
                    bottom=True, left=True, right=True, top=True)

c = ax.contourf(x,y,precip_mean.T,
                levels=np.arange(200,3300,100),
                extend='both',
        transform=ccrs.PlateCarree(),
          cmap=cmaps.NCV_jet)
cb=plt.colorbar(c,
                shrink=0.98,
                orientation='vertical',
                aspect=28,
                )
cb.ax.tick_params(labelsize=10,which='both',direction='in',)
plt.show()
fig.savefig('./TRMM_10year_monthly.png',dpi=500)

全部代码

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri May  5 09:45:08 2023

@author: %(jixianpu)s

Email : 211311040008@hhu.edu.cn

introduction : keep learning althongh walk slowly
"""
import xarray as xr
import os,glob
import numpy as np
from pyhdf.SD import SD
path = '/Datadisk/TRMM/3B43/'
file_list = []
for year in range(2009,2019):
    
    folder = os.path.join(path, str(year))
    
    file_name = glob.glob(folder+'/3B43.'+str(year)+'*.HDF')
    file_name.sort()
    file_list.extend(file_name)
    
file_list.sort()
dt = xr.open_dataset("/gpcp_monthly_mask.nc")
precip = dt.sel(time=slice('2009','2018'),lat =slice(-10,55),lon = slice(90,145)).precip
time_num = precip.time.values

def get_data(path,time):
    da = SD(path)
    pre = da.select('precipitation')[:]
    pre = np.expand_dims(pre,axis=0)
    lon = np.arange(-180,180.,0.25)
    lat = np.arange(-50,50.,0.25)
    time = time
#   time = datetime.datetime.utcfromtimestamp(tim).strftime('%Y-%m-%d %H:%M:%S')
    da = xr.DataArray(pre, 
                          dims=['time','lon','lat'],
                          coords=dict(
                              lon=(['lon'], lon),
                              lat=(['lat'], lat),
                              time=(['time'],[time])),
                              )
##############################################################################
    lon_range = [90.0, 145]
    lat_range = [-10, 55]
    da = da.sel(lon=slice(*lon_range), lat=slice(*lat_range))
    x ,y = da.lon,da.lat    
    return da,x,y

rain = np.zeros((len(file_list),221,240))          
for i in range(len(file_list)):
    print(i)
    da,x,y =  get_data(file_list[i], time_num[i])
    rain[i] = da
ds = xr.DataArray(rain, 
                      dims=['time','lon','lat'],
                      coords=dict(
                          lon=(['lon'], x.data),
                          lat=(['lat'], y.data),
                          time=(['time'],time_num)),
                          )
from global_land_mask import globe
def mask_land(ds, label='land', lonname='lon'):
    if lonname == 'lon':
        lat = ds.lat.data
        lon = ds.lon.data
        if np.any(lon > 180):
            lon = lon - 180
            lons, lats = np.meshgrid(lon, lat)
            mask = globe.is_ocean(lats, lons)
            temp = []
            temp = mask[:, 0:(len(lon) // 2)].copy()
            mask[:, 0:(len(lon) // 2)] = mask[:, (len(lon) // 2):]
            mask[:, (len(lon) // 2):] = temp
        else:
            lons, lats = np.meshgrid(lon, lat)# Make a grid
            mask = globe.is_ocean(lats, lons)# Get whether the points are on ocean.
        ds.coords['mask'] = (('lat', 'lon'), mask)
    elif lonname == 'longitude':
        lat = ds.latitude.data
        lon = ds.longitude.data
        if np.any(lon > 180):
            lon = lon - 180
            lons, lats = np.meshgrid(lon, lat)
            mask = globe.is_ocean(lats, lons)
            temp = []
            temp = mask[:, 0:(len(lon) // 2)].copy()
            mask[:, 0:(len(lon) // 2)] = mask[:, (len(lon) // 2):]
            mask[:, (len(lon) // 2):] = temp
        else:
            lons, lats = np.meshgrid(lon, lat)
            mask = globe.is_ocean(lats, lons)
        lons, lats = np.meshgrid(lon, lat)
        mask = globe.is_ocean(lats, lons)
        ds.coords['mask'] = (('latitude', 'longitude'), mask)
    if label == 'land':
        ds = ds.where(ds.mask == True)
    elif label == 'ocean':
        ds = ds.where(ds.mask == False)
    return ds
data =  mask_land(ds,'land')
precip_mean = np.nanmean(data*24*365,axis=0)
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import cmaps
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cnmaps import get_adm_maps, draw_maps
box = [100,140,0,50]
xstep,ystep = 10,10
proj = ccrs.PlateCarree(central_longitude=180)
plt.rcParams['font.family'] = 'Times New Roman',
fig = plt.figure(figsize=(8,7),dpi=200)
fig.tight_layout()
ax = fig.add_axes([0.1,0.2,0.8,0.7],projection = proj)
ax.set_extent(box,crs=ccrs.PlateCarree())
draw_maps(get_adm_maps(level='国'))
# ax.coastlines('50m')
ax.set_xticks(np.arange(box[0],box[1]+xstep, xstep),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(box[2], box[3]+1, ystep),crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=False)#True/False
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_title('TRMM(mm/year)',fontsize=16,pad=8,loc='left')
ax.tick_params(    which='both',direction='in',
                width=0.7,
                    pad=8, 
                    labelsize=14,
                    bottom=True, left=True, right=True, top=True)

c = ax.contourf(x,y,precip_mean.T,
                levels=np.arange(200,3300,100),
                extend='both',
        transform=ccrs.PlateCarree(),
          cmap=cmaps.NCV_jet)
cb=plt.colorbar(c,
                shrink=0.98,
                orientation='vertical',
                aspect=28,
                )
cb.ax.tick_params(labelsize=10,which='both',direction='in',)
plt.show()
fig.savefig('./TRMM_10year_monthly.png',dpi=500)

引用参考

TRMM: Tropical Rainfall Measuring Mission
https://climatedataguide.ucar.edu/climate-data/trmm-tropical-rainfall-measuring-mission
在这里插入图片描述

Monthly 0.25° x 0.25° TRMM multi-satellite and Other Sources Rainfall (3B43)
http://apdrc.soest.hawaii.edu/datadoc/trmm_3b43.php
在这里插入图片描述

TRMM 3B43: Monthly Precipitation Estimates
https://developers.google.com/earth-engine/datasets/catalog/TRMM_3B43V7
在这里插入图片描述

中巴经济走廊TRMM_3B43月降水数据（1998-2017年）：http://www.ncdc.ac.cn/portal/metadata/4b9504fa-0e34-47c9-a755-91d3f3253312
在这里插入图片描述

TRMM (TMPA/3B43) Rainfall Estimate L3 1 month 0.25 degree x 0.25 degree V7 (TRMM_3B43)（GES 官网介绍）：https://disc.gsfc.nasa.gov/datasets/TRMM_3B43_7/summary
在这里插入图片描述

国家海洋遥感在线分析平台 https://www.satco2.com/index.php?m=content&c=index&a=show&catid=317&id=217
在这里插入图片描述

之前没怎么处理过程HDF的文件，时间仓促，只是简单了记录了一下，没有考虑代码的美观和计算的高效性，欢迎大家评论或者联系我进行交流讨论~~

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

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

刚转岗做项目经理，无从下手，怎么办？

刚转岗做项目经理，无从下手，怎么办？

01 背景 最近在知乎平台看到一个问题是这么说的： 或许很多人都不是从工作开始就是项目专员再到项目经理这里一步一步过来，而是从其他岗位比如售前、销售、产品经理、程序员等转到项目经理岗位的。 那么对于这些人来说，做项目经理会有什么问…
阅读更多...
第一次找实习, 什么项目可以给自己加分（笔记）

第一次找实习, 什么项目可以给自己加分（笔记）

什么样的项目能简历加分、对找工作有帮助 基本特征： 一个特征是“硬核基础软件”，另一个为很实用的APP。 硬核基础软件 独立实现一个操作系统的kerne内核（操作系统的内部引擎） 北美计算机名校会让学生用一个学期的时间实现一个…
阅读更多...
冒险岛私人服务器详细架设教程

冒险岛私人服务器详细架设教程

冒险岛Online 　　《冒险岛Online》是由韩国WIZET和NEXON制作开发的一款2D横版卷轴网络游戏，于2004年7月24日在中国大陆正式上线，由盛大游戏负责运营。 　　故事以被“黑暗力量”不断入侵，因而进入了“浑沌期”的世界为背景，勇士们…
阅读更多...
微服务知识3

微服务知识3

Gateway核心概念 路由(route) 网关中最基础的部分，路由信息包括一个ID，一个目的URI，一组谓词工厂，一组Filter组成。如果谓词为真，则说明请求的URL和配置的路由匹配 谓词(preducates) 即java.util.function.Predica…
阅读更多...
atbf中imu数据的读取与处理方式

atbf中imu数据的读取与处理方式

一、说明 本文为作者在阅读atbf源码的过程中，对atbf中imu数据的读取和处理方式的个人理解，可能存在不对之处，意在抛砖引玉，请各位老师多多指正； 二、数据读取流程图 1、target NEUTRONRCF435SE 不同的target所定义…
阅读更多...
C++中的queue与priority_queue

C++中的queue与priority_queue

文章目录 queuequeue的介绍queue的使用 priority_queuepriority_queue介绍priority_queue使用 queue queue的介绍 队列是一种容器适配器，专门用于上下文先进先出的操作中。队列的特性是先进先出，从容器的一端插入，另一端提取元素。   队列…
阅读更多...
Java17的性能优势是否足以让它取代Java8？

Java17的性能优势是否足以让它取代Java8？

随着时间的推移，Java不断地进行更新和发展，以满足不断变化的业务需求。目前，Java8已经成为了一个非常成熟的版本，并且在各个领域广泛应用。但是，Java17也早已发布，并且是Java 11以来又一个LTS(长期支持)版本…
阅读更多...
vmware虚拟机安装k8s（之前已经安装过docker）

vmware虚拟机安装k8s（之前已经安装过docker）

1、安装开始 先执行：curl https://mirrors.aliyun.com/kubernetes/apt/doc/apt-key.gpg | sudo apt-key add 再执行更改源：echo "deb https://mirrors.aliyun.com/kubernetes/apt kubernetes-xenial main" >> /etc/apt/sources.list …
阅读更多...
单向链表——C语言实现

单向链表——C语言实现

哈喽，大家好，今天我们学习的是数据结构里的链表，这里主要讲的是不带哨兵卫头节点的单向链表，下篇将会继续带大家学习双向链表。 目录 1.链表的概念 2.单向链表接口的实现 2.1动态申请一个节点 2.2单链表打印 2.3单链表尾插 …
阅读更多...
强大的editplus 5.7

强大的editplus 5.7

EditPlus是一款由韩国 Sangil Kim （ES-Computing）出品的小巧但是功能强大的可处理文本、HTML和程序语言的Windows编辑器，你甚至可以通过设置用户工具将其作为C,Java,Php等等语言的一个简单的IDE。 EditPlus（文字编辑器&#xff0…
阅读更多...
Java学习-MySQL-数据库的设计

Java学习-MySQL-数据库的设计

Java学习-MySQL-数据库的设计 为什么需要设计数据库 当数据库比较复杂的时候，需要设计数据库。 糟糕的数据库设计： 数据冗余，浪费空间；数据库插入和删除很麻烦，可能导致异常（屏蔽使用物理外键&#xff0…
阅读更多...
机器学习实战教程（十三）：集成学习

机器学习实战教程（十三）：集成学习

简介 集成学习是一种机器学习方法，它旨在通过将多个单独的学习器（称为基分类器或基学习器）的预测结果进行组合，来提高整体的预测准确率。 集成学习可以看作是一种“多个人一起合作做事”的方法。每个基分类器都是独立的学习器&a…
阅读更多...
【恭喜宿主：你的神装Xpath到手】——07全栈开发——如桃花来

【恭喜宿主：你的神装Xpath到手】——07全栈开发——如桃花来

目录索引 什么是XML：文档演示： XML的节点关系：1.父节点：2. 子节点：3. 同胞节点：4. 先辈节点：5. 后代节点： Xpath：1. 相关语法：*最常用的路径表达式&#xff1…
阅读更多...
Android Studio实现文艺阅读App

Android Studio实现文艺阅读App

项目目录 一、系统概述二、系统特点三、开发环境四、运行演示五、源码获取 一、系统概述 本次带来的文艺阅读App可以提供高质量的原创文学作品。用户可以App中找到各种类型的文学作品，包括小说、散文、诗歌等，由来自不同领域的作家所创作。此外&#xf…
阅读更多...
对不起，我们不招还在用Excel的人，和金山系新秀比起差太远了

对不起，我们不招还在用Excel的人，和金山系新秀比起差太远了

相信点进来的朋友曾经也深受Excel荼毒。 的确，现如今在网上随便一搜，关于Excel的学习资料和答疑解惑的帖子不胜枚举，盖因为Excel有时太过热心，当然，是帮倒忙的那种热心。 自动把天数转换为日期，还替你把身…
阅读更多...
Flume 从入门到精通

Flume 从入门到精通

Flume Flume 是一种分布式、可靠且可用的服务 高效收集、聚合和移动大量日志 数据。 它具有基于流媒体的简单灵活的架构 数据流。它坚固耐用，容错，可靠性可调 机制以及许多故障转移和恢复机制。 它 使用允许在线分析的简单可扩展数据模型 应用。 系统要求…
阅读更多...
Java程序设计入门教程---循环结构（while）

Java程序设计入门教程---循环结构（while）

目录 思考 概念 语法 案例：求1到100的整数和？ 案例分析 思考 1. 让你输出10000000000000000句“Hello,world!”，你怎么写代码？ 2. 求1到100的整数和？ 概念 循环结构程序多次循环执行相同或相近的任务。 while循环…
阅读更多...
Nacos注册中心一些配置说明

Nacos注册中心一些配置说明

未安装Nacos的可以参考以下的安装教程 Nacos 安装教程（史上最详细保姆级教程）_nacos安装_大三的土狗的博客-CSDN博客 注意: Nacos默认是集群部署,如果想单机启动需要在对应Nacos的bin目录执行下面的命令 因为以后不可能只有一个Nacos注册中心,所以就默认…
阅读更多...
BERT+TextCNN实现医疗意图识别项目

BERT+TextCNN实现医疗意图识别项目

BERTTextCNN实现医疗意图识别项目 一、说明 本项目采用医疗意图识别数据集CMID传送门 数据集示例： {"originalText": "间质性肺炎的症状?", "entities": [{"label_type": "疾病和诊断", "start_pos&quo…
阅读更多...
Qt-数据库开发-用户登录、后台管理用户

Qt-数据库开发-用户登录、后台管理用户

Qt-数据库开发-用户登录、后台管理用户 [1] Qt-数据库开发-用户登录、后台管理用户1、概述2、实现效果 [2] Qt使用SqlLite实现权限管理初始化数据库创建数据表插入数据可使用结构体对数据信息进行封装数据库查询函数为数据库更新数据函数为删除数据函数为 [3] 测试效果 [1] Qt-…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023