R语言读取大型NetCDF文件

失踪人口回归,本篇来介绍下R语言读取大型NetCDF文件的一些实践。

1 NetCDF数据简介

先给一段Wiki上关于NetCDF的定义。

NetCDF (Network Common Data Form) is a set of software libraries and self-describing, machine-independent data formats that support the creation, access, and sharing of array-oriented scientific data. The project homepage is hosted by the Unidata program at the University Corporation for Atmospheric Research (UCAR). They are also the chief source of netCDF software, standards development, updates, etc. The format is an open standard. NetCDF Classic and 64-bit Offset Format are an international standard of the Open Geospatial Consortium.The project started in 1988 and is still actively supported by UCAR. The original netCDF binary format (released in 1990, now known as "netCDF classic format") is still widely used across the world and continues to be fully supported in all netCDF releases. Version 4.0 (released in 2008) allowed the use of the HDF5 data file format. Version 4.1 (2010) added support for C and Fortran client access to specified subsets of remote data via OPeNDAP. Version 4.3.0 (2012) added a CMake build system for Windows builds. Version 4.7.0 (2019) added support for reading Amazon S3 objects. Further releases are planned to improve performance, add features, and fix bugs.

本质上NetCDF是一个多维矩阵的数据,常用于地球科学领域的数据存储。

wiki百科定义

给出一个典型的例子(CHAP的O3数据)。

2 R语言处理大型NetCDF文件

我们可以看到NetCDF本质上这就是多个栅格叠在一起的文件,在R里面的处理方式基本依赖于几个栅格和NetCDF相关的包。包括ncdf4, raster/terra。

install.packages('ncdf4')
install.packages('raster')
install.packages('terra')

接下来给出一个标准的读nc文件的代码。

#读取文件
nc_o3 <- nc_open('CHAP_O3_Y10K_2013_V1.nc')#显示文件内容
print(nc_o3)#获取经纬度,变量名与缺失值
long <- ncvar_get(nc_o3, 'lon')
lat <- ncvar_get(nc_o3, 'lat')
o3 <- ncvar_get(nc_o3, 'O3')
fillvalue <- ncatt_get(nc_o3, "O3", "fill_value")
o3[o3==fillvalue$value] <- NA#生成栅格
r <- raster(t(o3), xmn=min(long), xmx=max(long), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +datum=WGS84 +no_defs"))#绘图
plot(r)

这部分代码的运行结果就是第一部分显示的那张图。这里要注意的是,print(nc_o3)那句代码是会决定下面获取经纬度,变量名与缺失名的关键。如图所示,这里的变量数量是一个,就是臭氧浓度O3,单位是 μ g / m 3 \mu g/m^3 μg/m3,变量名叫o3,缺失值为-999,然后对应的经纬度名字是lat和lon。

由于这个数据是10 km的处理起来比较快。当遇到全球或者全国比较精细化的NetCDF文件的时候,读取和另存为栅格可能会非常耗内存,因为R语言在处理数据的时候,默认是把数据全部读进去内存。笔者最近处理了一个全国km级的逐月气象数据。当我加载某一个三年的数据时,忽然内存飙升了,存进去的多维矩阵能占据5.7G内存。因此这就对处理速度和机子要求很高,也会有很多麻烦。那么当然我并不需要全国的数据,实际上我也是裁出研究区的数据。因此做了一下搜索和实践,总结了下如何根据需求,只读取部分区域的大型NetCDF文件。这样子就不需要机子内存的高要求了,这里以福建省为例(福建省的shp数据可以参照如下的链接下载)。使用的NetCDF数据为国家青藏高原科学数据中心1901到2022年逐月1km降水数据。

福建省标准地图矢量数据首次公开发布

中国1km分辨率逐月降水量数据集(1901-2022)

中国1km分辨率逐月平均气温数据集(1901-2022)

主要的代码如下:

#install.packages('sf')
#install.packages('raster')
#install.packages('terra')
#install.packages('ncdf4')
#install.packages('rasterVis')
#Load library
library(sf)
library(ncdf4)
library(raster)
library(terra)
library(rasterVis)#Read shapefile
fjcities <- read_sf('city.shp')
fjcitieswgs <- st_transform(fjcities, crs = '+proj=longlat +datum=WGS84 +no_defs')#Read netcdf file
nc_pre <- nc_open('pre_2021.nc')#Display the information of NetCDF
print(nc_pre)fjcitieswgs#To obtain the longitude, latitude, time, and name of variable
long <- ncvar_get(nc_pre, 'lon')
lat <- ncvar_get(nc_pre, 'lat')#Calculate numbers of rows and columns of the specific study area
LonIdx <- which(long[] > 115 & long[] < 121)
LatIdx <- which(lat[] > 23 & lat[] < 29)pre <- ncvar_get(nc_pre, 'pre', start = c(LonIdx[1], LatIdx[1],1),count = c(length(LonIdx), length(LatIdx),1))
fillvalue <- ncatt_get(nc_pre, "pre", "missing_value")
pre[pre==fillvalue$value] <- NA#To generate raster
r <- raster(t(pre), xmn=min(long[LonIdx]), xmx=max(long[LonIdx]), ymn=min(lat[LatIdx]), ymx=max(lat[LatIdx]), crs=CRS("+proj=longlat +datum=WGS84 +no_defs"))#Display the raster
plot(r)#Clip the raster using the shapefile
r <- crop(r, fjcitieswgs)
r <- mask(r, fjcitieswgs)#Display the raster
plot(r)#Using different r package to plot raster
levelplot(r)

关键在于ncvar_get函数里的start和count参数。这两个负责控制读取NC的行列数据以及多维数(如果没有时间轴,直接给一个2个元素的数组就行)。start是读取NetCDF的起始行列数,count是读取NetCDF的数量。后续转栅格的操作是raster函数的写法。由于raster快停止维护了。我们也会提供terra包的写法(实际上差别不大)。

#Use terra pacakge to generate raster
rn <- rast(t(pre), crs= '+proj=longlat +datum=WGS84 +no_defs')#Display the raster
plot(rn)

由于后续的裁剪代码terra和raster毫无差别。这里就不赘述了。这部分代码我也会放在我的Github项目上(My Studies of Urban GIS)。

最后感谢下Google搜到的几个资料。

参考链接:

  • Read multiple layers raster from ncdf file using terra package

  • extracting large layer netcdf using r

  • How to reduce the size of a NetCDF in R?

  • Extracting a large layer from NetCDF using R

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

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

GlobalExceptionHandler全局异常处理器的设计

在Java Web开发中&#xff0c;全局异常处理器&#xff08;GlobalExceptionHandler&#xff09;是一个非常重要的概念。它允许我们集中处理应用程序中可能发生的各种异常&#xff0c;从而提供统一的错误响应&#xff0c;增强用户体验&#xff0c;并简化异常处理逻辑。下面将详细…

栈队列数组试题(二)——队列

一、单项选择题 01.栈和队列的主要区别在于&#xff08;). A.它们的逻辑结构不一样 B.它们的存储结构不一样 C.所包含的元素不一样 D.插入、删除操作的限定不一样 02&#xff0e;队列的“先进先出…

STM32串口通信—串口的接收和发送详解

目录 前言&#xff1a; STM32串口通信基础知识&#xff1a; 1&#xff0c;STM32里的串口通信 2&#xff0c;串口的发送和接收 串口发送&#xff1a; 串口接收&#xff1a; 串口在STM32中的配置&#xff1a; 1. RCC开启USART、串口TX/RX所对应的GPIO口 2. 初始化GPIO口 …

Linux mktemp命令教程:创建临时文件和目录(附实例详解和注意事项)

Linux mktemp命令介绍 mktemp命令在Linux中用于创建临时文件或目录。这个命令可以帮助我们在脚本或应用程序中创建一个有效且唯一的临时文件或目录。 Linux mktemp命令适用的Linux版本 mktemp命令在所有主要的Linux发行版中都可以使用&#xff0c;包括Debian、Ubuntu、Alpin…

鸿蒙跨包跳转页面-HSP页面路由

页面路由跳转 若开发者想在entry模块中&#xff0c;添加一个按钮跳转至library模块中的menu页面&#xff08;路径为&#xff1a;library/src/main/ets/pages/menu.ets&#xff09;&#xff0c;那么可以在使用方的代码&#xff08;entry模块下的Index.ets&#xff0c;路径为&am…

YOLOv8改进 | 图像去雾 | 特征融合注意网络FFA-Net增强YOLOv8对于模糊图片检测能力(北大和北航联合提出)

一、本文介绍 本文给大家带来的改进机制是由北大和北航联合提出的FFA-net: Feature Fusion Attention Network for Single Image Dehazing图像增强去雾网络,该网络的主要思想是利用特征融合注意力网络(Feature Fusion Attention Network)直接恢复无雾图像,FFA-Net通过特征…

MyBatis-Plus学习记录

目录 MyBatis-Plus快速入门 简介 快速入门 MyBatis-Plus核心功能 基于Mapper接口 CRUD 对比mybatis和mybatis-plus&#xff1a; CRUD方法介绍&#xff1a; 基于Service接口 CRUD 对比Mapper接口CRUD区别&#xff1a; 为什么要加强service层&#xff1a; 使用方式 CR…

开发指南009-从list导出excel文件

从数据库返回一般是对象的列表&#xff0c;平台底层提供了从list转为excel文件的方法。平台的设计思想就是为一些典型的场景设计对应的解决方法&#xff0c;通过模式化的方法来简化编程和提高维护性&#xff08;通过标准化来减少学习成本和维护成本&#xff0c;张三做的东西和李…

ARC 的 retainCount 是怎么存储的

ARC的retainCount是存吃在64张哈希表中的&#xff0c;根据哈希算法去查找所在的位置&#xff0c;无需便利 散列表&#xff08;引用计数表、weak表&#xff09; SideTables 表在 非嵌入式的64位系统中&#xff0c;有 64张 SideTable 表每一张 SideTable 主要是由三部分组成。自…

【MASM汇编语言快速入门】MASM常用伪指令速查表——变量

MASM伪指令速查表–变量 初学MASM时&#xff0c; 常常看不懂db, dup(?)等汇编指令的含义&#xff0c; 教材中也缺乏系统的解释。与机器指令不同&#xff0c;这些指令叫伪指令&#xff0c; 在编译&#xff08;汇编&#xff09;的时候被MASM编译器处理&#xff0c; 而在运行时计…

Boot——组件(导航和选项卡、分页、卡片、轮播图)

Boot——组件&#xff08;下&#xff09; 导航和选项卡 https://v5.bootcss.com/docs/components/navs-tabs/ &#xff08;1&#xff09;导航 <ul class"nav"><li class"nav-item"><a href"#" class"nav-link">…

游戏数据处理

游戏行业关键数据指标 ~ 总激活码发放量、总激活量、总登录账号数 激活率、激活登录率 激活率 激活量 / 安装量 激活率 激活量 / 激活码发放量 激活且登录率 激活且登录量 / 激活码激活量 激活且登录率应用场景 激活且登录率是非常常用的转化率指标之一&#xff0c;广泛…

Ypay源支付6.9无授权聚合免签系统可运营源码

YPay是一款专为个人站长设计的聚合免签系统&#xff0c;YPay基于高性能的ThinkPHP 6.1.2 Layui PearAdmin架构&#xff0c;提供了实时监控和管理的功能&#xff0c;让您随时随地掌握系统运营情况。 说明 Ypay源支付6.9无授权聚合免签系统可运营源码 已搭建测试无加密版本…

HTML5:七天学会基础动画网页13

看完前面很多人可能还不是很明白0%-100%那到底是怎么回事&#xff0c;到底该怎么用&#xff0c;这里我们做一个普遍的练习——心跳动画 想让心❤跳起来&#xff0c;我们先分析一波&#xff0c;这个心怎么写&#xff0c;我们先写一个正方形&#xff0c;再令一个圆形前移: 再来一…

蓝桥杯历年真题 省赛 Java b组 2016年第七届

一、题目 分小组 9名运动员参加比赛&#xff0c;需要分3组进行预赛。 有哪些分组的方案呢&#xff1f; 我们标记运动员为 A,B,C,... I 下面的程序列出了所有的分组方法。 该程序的正常输出为&#xff1a; ABC DEF GHI ABC DEG FHI ABC DEH FGI ABC DEI FGH ABC DFG EHI ABC…

Linux中YUM仓库的配置

Linux软件包的管理 YUM仓库是什么YUM的常用命令修改YUM源其实CentOS7已经对YUM做了优化 YUM仓库是什么 之前传统RPM的管理方式 可以简单理解为写Java的时候不用Maven管理 jar包都要自己手动去导入 去下载 但是配置好YUM仓库 就放佛在用Maven管理Java项目 基于RPM包管理 能够从…

Python导入类说一说

要在Python中导入一个类&#xff0c;需要使用import关键字。 详细去看下面的代码 1、多例类 class Restaurant:餐馆类def __init__(self,restaurant_name,cuisine_type):#类的属性self.restaurant_name restaurant_nameself.cuisine_type cuisine_type# self.stregth_leve…

2024软件测试应该学什么?“我“怎么从功能转入自动化测试?

目录&#xff1a;导读 前言一、Python编程入门到精通二、接口自动化项目实战三、Web自动化项目实战四、App自动化项目实战五、一线大厂简历六、测试开发DevOps体系七、常用自动化测试工具八、JMeter性能测试九、总结&#xff08;尾部小惊喜&#xff09; 前言 1、软件测试应该学…

TypeScript编译选项

编译单个文件&#xff1a;终端 tsc 文件名 自动编译单个文件&#xff1a;终端 tsc 文件名 -w 编译整个项目&#xff1a;tsc 前提是得有ts的配置文件tsconfig.json 自动编译整个项目&#xff1a;tsc --w tsconfig.json默认文件内容&#xff1a; tsconfig.json是ts编译器的配…

代码随想录算法训练营第16天

104.二叉树的最大深度 &#xff08;优先掌握递归&#xff09; 思路&#xff1a; 注意&#xff1a; 传入参数&#xff1a;depth, root 终止条件&#xff1a;if(root nullptr) return 0; 单层递归逻辑&#xff1a; 左右中int left getmax(depth1, root->left);int right …