我在读取一个nc数据的时候,想通过读取一系列用户自定义的经纬度坐标来获取计算出来的温度数据 。但是,我遇到一个问题:nc数据本身没有经纬度信息!数据本身的情况如下:
nc数据储存的温度数据TBOT信息有三维:time=124, lat=360, lon=720。因为储存了每天每6小时的温度数据,一天有4个6小时,所以一个月一共有31*4=124,即time=124。 lat代表纬度,即从南纬-90°到北纬90°,每0.5°为一个步长储存温度数据,所以lat=(90+90)/0.5=360。 lon代表经度,即从0°到360°,每0.5°为一个步长储存温度数据,所以lat=360/0.5=720。
因为三维数据打不开,我对数据进行了重塑压缩,形成259200x124 double的矩阵t_2d,并进行了后续一系列计算。259200=720*360,相当于用一系列经纬度索引确定了一个位置,这个大家可以理解吧。
问题来了!259200这每一行代表的是地球上每一点的位置的温度值,但是这个位置是从哪里开始,朝着什么方向,在哪里结束呢?我猜测是从(0.25 ,-89.75)这个地方开始到(359.75,89.75)结束。所以,我想通过计算输入的经纬度相对于起始位置(0.25 ,-89.75)的偏移量来找到指定经纬度值对应的索引,然后将这个索引的经纬度相乘,得到的最终数值就是该位置对应259200中的行数。可是我发现这似乎不对。我发现一些地方的温度与其季节明显不符,一看就知道错的。所以这个问题,应该怎么解决
clc
clear all;
% 设置循环参数
start_month = 7;
end_month = 8;
% 存储每日和每月各温度指标的变量
avg_t_daily = [];
max_temp_daily= [];
min_temp_daily= [];
avg_t_monthly = [];
max_t_monthly = [];
min_t_monthly = [];
temperature=[];
studyaera=[];
LA=[];
% 循环读取和处理多个nc图像
for month = start_month:end_month
% 构建nc图像文件名
file_name = sprintf('F:\\TPHWL6Hrly\\clmforc.cruncep.V7.c2016.0.5d.TPQWL.2015-%02d.nc', month);
% 读取温度数据
t = ncread(file_name, 'TBOT');
LONGXY= ncread(file_name, 'LONGXY');
LATIXY= ncread(file_name, 'LATIXY');
% 获取变量 t 的尺寸信息(依次为:经度、纬度、时间),并将其分别存储在自定义变量 num_lon、num_lat 和 num_time
[num_lon, num_lat, num_time] = size(t);
% 重塑矩阵为二维
t_2d = reshape(t, [], num_time);
%% 求每天的气温情况,包括平均气温,日高温,日低温,365个值
% 计算365天每天的平均气温
num_days = floor(num_time / 4); % 总共的天数,向下取整
avg_t_daily_month = zeros(size(t_2d, 1), num_days); % 存储每天的平均温度
for i = 1:num_days
start_index = (i - 1) * 4 + 1;
end_index = i * 4;
avg_t_daily_month(:, i) = mean(t_2d(:, start_index:end_index), 2) - 273.15;
end
% 将每日平均温度存储到整体变量中
avg_t_daily = [avg_t_daily avg_t_daily_month];
% 计算365天每天的最高温度
max_temp_daily_month = zeros(size(t_2d, 1), num_days); % 存储每天的最高温度
for i = 1:num_days
start_index = (i - 1) * 4 + 1;
end_index = i * 4;
max_temp_daily_month(:, i) = max(t_2d(:, start_index:end_index), [], 2) - 273.15;
end
% 将365天每天最高温度存储到整体变量中
max_temp_daily = [max_temp_daily max_temp_daily_month];
% 计算365天每天的最低温度
min_temp_daily_month = zeros(size(t_2d, 1), num_days); % 存储每天的最高温度
for i = 1:num_days
start_index = (i - 1) * 4 + 1;
end_index = i * 4;
min_temp_daily_month(:, i) = min(t_2d(:, start_index:end_index), [], 2) - 273.15;
end
% 将计算365天每天的最低温度存储到整体变量中
min_temp_daily = [min_temp_daily min_temp_daily_month];
end
%% 求每月的气温情况,包括平均气温,月高温,月低温,12个值
% 计算每月平均气温
avg_t_mons_month = mean(t_2d, 2) - 273.15; % 沿纬度方向求平均值
% 将每个月的平均温度存储到整体变量中
avg_t_monthly = [avg_t_monthly avg_t_mons_month];
% 计算每月最高温
max_t_mons_month = max(t_2d, [], 2) - 273.15; % 沿纬度方向求平均值
% 将每个月的最高温存储到整体变量中
max_t_monthly = [max_t_monthly max_t_mons_month];
% 计算每月最低温
min_t_mons_month = min(t_2d, [], 2) - 273.15; % 沿纬度方向求平均值
% 将每个月的最低温存储到整体变量中
min_t_monthly = [min_t_monthly min_t_mons_month];
%% 求自定义地区的气温情况
longitudes = [113.17, 47.49, 110.75];
latitudes = [23.8, 30.30, 22.75];
% 定义纬度和经度的范围
lat_range = [-89.75, 89.75];
lon_range = [0.25, 359.75];
% 计算纬度和经度的步长
lat_step = 0.5;
lon_step = 0.5;
% 找到指定经纬度值对应的索引
lat_indices = floor((latitudes - lat_range(1)) / lat_step) + 1;
lon_indices = floor((longitudes - lon_range(1)) / lon_step) + 1;
% 经纬度值组合最终索引
result = lat_indices + (lon_indices - 1) * (lat_range(2) - lat_range(1) + 1) / lat_step;
% 根据最终索引循环获取对应的 t_2d 值
t_2d_studyarea = t_2d(result, :)-273.15;
avg_t_daily_studyarea = avg_t_daily(result, :)
min_temp_daily_studyarea = min_temp_daily(result, :)
max_temp_daily_studyarea = max_temp_daily(result, :)
avg_t_mons_month_studyarea = avg_t_mons_month(result, :)
max_t_mons_month_studyarea = max_t_mons_month(result, :)
min_t_mons_month_studyarea = min_t_mons_month(result, :)
这要看你的 nc 数据是如何组织的,,我猜你可能是把这个混淆了,
如果数据是按纬度(lat)先变化,然后是经度(lon)变化的顺序存储的(就是每个纬度下存储了所有经度的数据),那么索引计算应该是:
result = (lon_indices - 1) * num_lat + lat_indices;
反之
result = (lat_indices - 1) * num_lon + lon_indices;
ersal Transverse Mercator)是一种平面直角坐标系,常用于地图制图和大地测量。在处理地理信息数据时,常需要将经纬度坐标转换为UTM坐标。
对于该问题,可以使用经纬度数据进行反向计算,得到每个位置对应的行号。具体解决方案如下:
lat_start = -89.75; % 起始纬度 lon_start = 0.25; % 起始经度 lat_delta = 0.5; % 纬度步长 lon_delta = 0.5; % 经度步长
lat_offset = floor((lat - lat_start) / lat_delta) + 1; lon_offset = floor((lon - lon_start) / lon_delta) + 1;
需要注意的是,因为数据中的纬度是从-90°到90°,而实际的纬度范围是从-90°到+90°,所以需要将计算得到的行索引加上一。
index = (lat_offset - 1) * 720 + lon_offset;
最终得到的index即为该位置在矩阵中的行号。这个行号是从1开始的,所以需要将其减去1再作为索引来获取相应位置的温度值。
下面是完整的MATLAB代码示例:
%% 定义经纬度信息 lat_start = -89.75; % 起始纬度 lon_start = 0.25; % 起始经度 lat_delta = 0.5; % 纬度步长 lon_delta = 0.5; % 经度步长
%% 读取nc文件的温度数据矩阵,假设为T T = rand(259200, 124); % 这里用随机数代替真实数据
%% 定义需要查询的经纬度坐标,假设为(40.5, 120) lat = 40.5; lon = 120;
%% 计算该位置相对于起始位置的偏移量 lat_offset = floor((lat - lat_start) / lat_delta) + 1; lon_offset = floor((lon - lon_start) / lon_delta) + 1;
%% 计算该位置在矩阵中的行号 index = (lat_offset - 1) * 720 + lon_offset - 1;
%% 获取该位置的温度值 temp = T(index, :); % 获取该位置在所有时间点上的温度值