将xyz格式的GRACE数据转成geotiff格式

我们需要将xyz格式的文件转成geotiff便于成图,或者geotiff转成xyz用于数据运算,下面介绍如何实现这一操作,采用GMT和matlab两种方法。

1.GMT转换

我们先准备一个xyz文件,这里是一个降水文件。在gmt中采用以下的语句实现xyz转grd网格文件

    xyz2grd DJF.txt -R-180/180/-90/90 -JN12c -I5 -Gm1.grdgmt grdsample m1.grd -Gm11.grd -I0.5

语句实现了将xyz转成grd文件,然后在Global mapper中打开,导出为geotiff格式即可,然后可以用于绘图和地理信息处理显示。

2.matlab程序实现

我们现在matlab中读取一个geotiff文件以观察其组成要素。

[A1,R] = geotiffread('Tibet_1000.tif');A1(A1<=0)=0;

可以看到,其由两部分组成,一个是数值矩阵(A1),另一个是投影信息(R)。因此我们如果想要将xyz数据转成geotiff格式的文件,同样需要准备两个信息:一个是数据矩阵,另一个是投影信息(当然需要自己设置)。

[A1,R] = geotiffread('Tibet_1000.tif');A1(A1<=0)=0;
lon = linspace(65.0042,109.9958,5400);
lat = linspace(20.0042,49.9958,3600);
[lon,lat] = meshgrid(lon,lat);O1.lon = lon;O1.lat = lat;O1.rg  = flipud(double(A1));
rg_plot(O1),title('10 km DEM'),colorbar

需要注意的是,一般我们所说的xyz数据是三列的,分别是

如果要转成geotiff文件,需要提前转成数值矩阵,即通常需要reshape一下(这针对于转换的影响为矩形)

接下来,我们实现一个xyz转成geotiff的例子。

(1)先准备一个xyz数据,这里以MSSA插值的GRACE level03数据集为例,数据参考以下专栏:

分享一个月份连续的MSSA插值的GRACE level03数据集_grace数据集-CSDN博客

matlab先读取xyz数据,然后绘图,可见reshape的正确性:

% % load data
A = load('GRACE_MSSA_2022_01.xyz');%% reshape
O.lon = reshape(A(:,1),181,361);
O.lat = reshape(A(:,2),181,361);
O.rg  = reshape(A(:,3),181,361);wzq_plot(O)

wzq_plot函数如下,其中缺失的报错文件参加B站的置顶评论:

绘图函数的使用wzq_plot - 哔哩哔哩 (bilibili.com)

function wzq_plot(wzq)
pcolor(wzq.lon,wzq.lat,wzq.rg)
shading interp
hold on;
if max(wzq.lon(:)) < 200coast=load('coastline-from-GMT-WNI.dat');
elsecoast=load('coastline-from-GMT-WNI-0-360.dat');
end
plot(coast(:,1),coast(:,2),'k')
hold off;
colorbar
colormap jet
end 

(2)配置投影信息

这里我们借用其他读取的geotiff文件的投影信息,然后按照实际情况进行修改配置,这里我们采用一个DEM的投影信息,我们需要修改的地方包括:

X Y的范围  采样分辨率  经纬度范围  栅格数量

然后我们得到配置好的新的投影信息,但是实际上有更简单的配置方法:

[A1,R] = geotiffread('Tibet_1000.tif');A1(A1<=0)=0;
% % load data
A = load('GRACE_MSSA_2022_01.xyz');%% reshape
O.lon = reshape(A(:,1),181,361);
O.lat = reshape(A(:,2),181,361);
O.rg  = reshape(A(:,3),181,361);wzq_plot(O)%% setting new projection info
R1.RasterInterpretation      = 'Postings';
R1.XIntrinsicLimits          = [1,361];
R1.YIntrinsicLimits          = [1,181];
R1.SampleSpacingInLatitude   = 1;
R1.SampleSpacingInLongitude  = 1;
R1.LatitudeLimits            = [-90,90];
R1.LongitudeLimits           = [0,360];
R1.RasterSize                = [181,361];
R1.AngleUnit                 = 'degree';
R1.ColumnsStartFrom          = 'north';
R1.RowsStartFrom             = 'east';
R1.CoordinateSystemType      = 'geographic';
R1.AngleUnits                = 'degrees';latlim = [-90,90];
lonlim = [0,360];
R0 = georefcells(latlim,lonlim,size(O.rg));	% 设置一个地理坐标
geotiffwrite('GRACE_xyz2tiff.tif', O.rg, R0);		% 

然后我们可以得到GRACE_xyz2tiff.tif文件,需要注意的是,再次运行前需要删除之前生成的文件。

在global mapper中可以打开

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

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

相关文章

智慧物业 内外共建 ——七巧低代码打造终端智能物业管理系统

项目背景 物业管理是城市管理的重要组成部分&#xff0c;也是社区居民的刚需服务。 随着城市化进程的加快&#xff0c;物业管理面临着越来越多的挑战&#xff0c;如业主需求的多样化、业务流程的复杂化、数据信息的庞大化、外部系统的集成化等。这些挑战要求物业公司提高自身…

一年发3篇TRO的“刺头”是怎样炼成的?

“世界上最成功的机器人,是人类本身。” 香港城市大学MetaSlam与GAIRLAB的创始人殷鹏教授,在知乎一则「世界上最厉害的机器人是什么」的问题下,写下了这样一句回答。 从小就痴迷于机器人的殷鹏,在2023这一年的时间里带领团队发表了3篇T-RO,这位90后的年轻教授展现出他在…

STM32定时器中断

定时器可以对输入的时钟进行计数&#xff0c;并在计数值达到设定值时发出中断 定时器就是一个计数器 预分频器&#xff1a;对系统时钟进行分频得到定时器时钟频率 自动重装在值&#xff1a;计数多少个进入中断 基本定时器两个&#xff0c;tim6和7&#xff0c;挂载在apb1 通…

【开源】基于JAVA+Vue+SpringBoot的软件学院思政案例库系统

目录 一、摘要1.1 项目介绍1.2 项目录屏 二、功能模块2.1 系统管理员2.2 普通教师 三、系统展示四、核心代码4.1 查询思政案例4.2 审核思政案例4.3 查询思政课程4.4 思政案例点赞4.5 新增思政案例评语 五、免责说明 一、摘要 1.1 项目介绍 基于JAVAVueSpringBootMySQL的软件学…

Linux---信号

前言 到饭点了&#xff0c;我点了一份外卖&#xff0c;然后又开了一把网游&#xff0c;这个时候&#xff0c;我在打游戏的过程中&#xff0c;我始终记得外卖小哥会随时给我打电话&#xff0c;通知我我去取外卖&#xff0c;这个时候游戏还没有结束。我在打游戏的过程中需要把外…

从[redis:LinkedList]中学习链表

文章目录 adlistlistNodelistmacros[宏定义]listCreatelistInitNodelistEmptylistReleaselistAddNodeHeadlistLinkNodeHeadlistAddNodeTaillistLinkNodeTaillistInsertNodelistDelNodelistUlinkNodelistIndexredis3.2.100quicklistredis7.2.2quicklist redis的基本数据类型之一…

day 19 (进阶)

一 首先 昨日内容回顾 思维导图&#xff1a;&#xff08;日更附 养成习惯 加油&#xff09; 补充Linux思维导图 衔接一下之前学过的 二 课堂知识提炼 练习&#xff1a;统计文件行数 想查看是否正确就用 grep -c “文件名” 来看 会输出结果 练习&#xff1a;把file.c里面的…

JS第二天、原型、原型链、正则

☆☆☆☆ 什么是原型&#xff1f; 构造函数的prototype 就是原型 专门保存所有子对象共有属性和方法的对象一个对象的原型就是它的构造函数的prototype属性的值。prototype是哪来的&#xff1f;所有的函数都有一个prototype属性当函数被创建的时候&#xff0c;prototype属性…

从资深用户角度谈三款出色数据可视化工具

作为一名数据可视化领域的老用户&#xff0c;我接触过众多数据可视化产品&#xff0c;其中不乏佼佼者。今天&#xff0c;我想为大家介绍三款在我心目中颇具特色的数据可视化产品&#xff0c;它们分别是山海鲸可视化、Tableau和Power BI。 首先&#xff0c;让我们来谈谈山海鲸可…

Python字符串和日期时间格式转换

Python字符串和日期时间格式转换 前言&#xff1a;1.字符串和日期时间转换终论&#xff1a;给定月份的上月月份(YYYYMM)1.1格式YYYYMM变成YYYYMMDD1.2字符串转换为时间格式1.3时间格式加减1.4时间格式转换为字符串 2.Pandas的DataFrame时间格式转换 前言&#xff1a; 字符串转…

vulhub中Apache APISIX Dashboard API权限绕过导致RCE(CVE-2021-45232)

Apache APISIX是一个动态、实时、高性能API网关&#xff0c;而Apache APISIX Dashboard是一个配套的前端面板。 Apache APISIX Dashboard 2.10.1版本前存在两个API/apisix/admin/migrate/export和/apisix/admin/migrate/import&#xff0c;他们没有经过droplet框架的权限验证&…

精细管理药厂设备,制药机械设备管理平台系统助力生产提效

制药行业的复杂性要求对药品的品质和安全性进行严格控制&#xff0c;而这离不开高效管理各类机械设备。然而&#xff0c;随着制药企业规模的不断扩大和技术的迅猛进步&#xff0c;如何有效管理这些设备成为一个亟待解决的问题。在这一挑战面前&#xff0c;PreMaint制药机械设备…