weighted variable on a satellite
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Multiplying the pixel values by the cosine of the latitude helps to account for this change in grid cell area. This scaling is often referred to as "area weighting" or "cosine weighting." It ensures that the contributions of grid cells at different latitudes are appropriately weighted in calculations of regional or global averages, taking into account the variation in grid cell size.
how can i do this to this code (i edited this code so that its matches the original)
close all;
clc;
clear all
format long g;
input = 'C:\Users\vasco\OneDrive\Ambiente de Trabalho\ESTAGIO_Vasco\JPL TELLUS GRACE Level-3\';
output = 'C:\Users\vasco\OneDrive\Ambiente de Trabalho\ESTAGIO_Vasco\JPL TELLUS GRACE Level-3 Output\';
lista_ficheiros = dir(fullfile(input, '*.nc'));
lon_1 = double(ncread(fullfile(input, lista_ficheiros(1).name), 'lon'));
lat_1 = double(ncread(fullfile(input, lista_ficheiros(1).name), 'lat'));
nlon = length(lon_1);
nlat = length(lat_1);
ntime = length(lista_ficheiros);
lwe_data = zeros(nlon, nlat, ntime);
for i = 1:length(lista_ficheiros)
filename = fullfile(input, lista_ficheiros(i).name);
nc = netcdf.open(filename);
lwe = ncread(filename, 'lwe_thickness');
lwe_data(:, :, i) = lwe;
netcdf.close(nc);
lat=ncread(filename, 'lat');
lon=ncread(filename, 'lon');
% longitude esta de [0->360]. transformo em [0->180, -180->0]
idx = lon>180;
lon(idx) = lon(idx) - 360;
% circshift lon para o primeiro elemento maior que 180 torna-se o primeiro elemento
% em lon tornando lon em [-180 -> 180]
avanco = 1-find(idx,1);
lon = circshift(lon,avanco);
% circshift lwe o mesmo valor na primeira dimensao sendo esta a dimensao lon
lwe = circshift(lwe,avanco,1);
time = double(ncread(filename, 'time'));
tm = regexprep(lista_ficheiros(i).name, '.*?(\d{7}-\d{7}).*', '$1');
nome_ficheiro = strcat('GRACE-', tm);
figure(1);
axis equal;
clev = -0.5:0.01:0.5;
[x,y] = meshgrid(lon,lat);
contourf(x, y, lwe', clev, 'LineStyle', 'none', 'Fill', 'on');
xlabel('Longitude');
ylabel('Latitude');
clim([min(clev), max(clev)]);
colormap(winter(length(clev)-1));
colorbar('eastoutside');
c = colorbar;
c.Label.String = 'Valor em metros (m)';
title(strcat('GRACE-', tm));
%estatisticas
nval = nnz(~isnan(lwe));
s = nansum(lwe);
media = nansum(s) / nval;
med_global(i)=media;
tempo(i)=time/365.25+2002.00;
%Guardar figuras
figura = fullfile(output, strcat(nome_ficheiro, '.png'));
saveas(gcf, figura);
estatisticas = fullfile(output, strcat(nome_ficheiro, '_estatisticas.txt'));
fid = fopen(estatisticas, 'w');
fprintf(fid, 'Estatisticas GRACE-%s:\n', tm);
fprintf(fid, 'Soma lwe_thickness: %f\n', s);
fprintf(fid, 'Media lwe_thickness: %f\n', media);
fclose(fid);
% nao fechei a figura pois da para ver a transicao de valores de uma
% para a outra
% close(gcf);
end
%%
%plot serie temporal - média regão imagem a imagem
figure(2);
title('Serie temporal')
%considerar anos inteiros
med_global2=med_global(1:158);
tempo2=tempo(1:158);
plot(tempo2,med_global2,'b--*')
grid on
xlabel('Tempo (anos)');
xticks(2002:1:2018);
ylabel('Média Global (m)');
nome_final_serie='Serie_temporal_GRACE';
hold on
%Ajuste linear
LWE_fit = polyfit(tempo2,med_global2,1);
slopeLWE=LWE_fit(1);
lin_ajusteLWE = polyval(LWE_fit,tempo2);
plot(tempo2,lin_ajusteLWE,'r', 'LineWidth', 1.5);
figura_2 = fullfile(output, strcat(nome_final_serie, '.png'));
saveas(gcf, figura_2);
%plot média da variável para o período todo, pixel a pixel
lwe_data_file = fullfile(output, 'lwe_data.mat');
save(lwe_data_file, 'lwe_data')
lwe_transpor = permute(lwe_data, [3, 1, 2]);
pixel_media = mean(lwe_transpor, 1);
pixel_std = std(lwe_data, 0, 3);
pixel_output_std = fullfile(output, 'pixel_std.mat');
pixel_outputs = fullfile(output, 'pixel_media.mat');
save(pixel_outputs, 'pixel_media');
nome_final='Media total da GRACE';
figure(3);
clevv = -0.1:0.01:0.1;
data = squeeze(pixel_media(1, :, :));
final = circshift(data,avanco,1);
contourf(lon, lat, final.', clevv, 'LineStyle', 'none', 'Fill', 'on');
xlabel('Longitude');
ylabel('Latitude');
title('Mapa da média de lwe');
clim([min(clevv), max(clevv)]);
c = colorbar;
c.Label.String = 'Valor em metros (m)';
colorbar('eastoutside');
colormap(hot(length(clevv) - 1));
hAx = gca;
hAx.YDir = 'normal';
hold on
C=load('coastlines');
plot(C.coastlon, C.coastlat, 'k')
figura_3 = fullfile(output, strcat(nome_final, '.png'));
saveas(gcf, figura_3);
nome_std='Mapa_std';
figure(4);
final2 = circshift(pixel_std,avanco,1);
contourf(lon, lat, final2.', clevv, 'LineStyle', 'none', 'Fill', 'on');
xlabel('Longitude');
ylabel('Latitude');
clevv = 0:0.01:0.2;
clim([min(clevv), max(clevv)]);
c = colorbar;
c.Label.String = 'Valor em metros (m)';
colorbar('eastoutside');
colormap(hot(length(clevv) - 1));
title=('Standard deviation');
hAx = gca;
hAx.YDir = 'normal';
hold on
C=load('coastlines');
plot(C.coastlon, C.coastlat, 'k')
figura_4 = fullfile(output, strcat(nome_std, '.png'));
saveas(gcf, figura_4);
0 comentarios
Respuestas (1)
Walter Roberson
el 10 de Sept. de 2023
clear all
we can be certain after that, that there are no "left-over" variables in the workspace -- that any variables not defined in the code are not expected to already exist.
for i = 1:length(lista_ficheiros)
There is no assignment to med_global before that loop.
The loop will be iterated once for each .nc file that was found in the appropriate directory.
med_global(i)=media;
Here med_global is assigned into. Because it was not initialized before the loop, med_global will grow in place, becoming only as large as the maximum number of iterations performed -- so med_global will have only one entry for each .nc file that was found in the appropriate directory.
% close(gcf);
end
There are no other assignments into med_global before the end of that loop.
med_global2=med_global(1:158);
Here, med_global is expected to be at least 158 entries long, implying that the code expects at least 158 .nc files to be present in the appropriate directory. But what if the directory only has one .nc file? What if the directory has more than 158 .nc files, but the 158 .nc files of interest do not happen to be the first 158 .nc files as returned by dir() .
Were you aware that dir() does not sort the file names in any way? dir() returns files in whatever order was present when it asked the operating system for information about the directory.
input = 'C:\Users\vasco\OneDrive\Ambiente de Trabalho\ESTAGIO_Vasco\JPL TELLUS GRACE Level-3\';
You are using Windows. Windows itself does not define any sort order for directories -- Windows uses whatever order the file system defines. Different file system types are completely permitted to use different orders -- including potentially "newest files always go at the end" or potentially keeping "slots" for files and writing new files into the first "unused" slot so the location for a new file might depend on the history of file transactions in the directory. If the file system involved turnes out to be NTFS, then NTFS does do sorting -- but the sort order for any given NTFS file system depends upon the Language setting of the person who created the NTFS file system (and does not depend upon the language setting of the user of the NTFS file system).
It is always a mistake to rely on the number and order of files returned by dir() without at least testing that the conditions hold. For example if you need at least 158 files, then immediately after you do the dir() test whether there at least 158 files present and give a clean error and do not execute the code.
10 comentarios
Ver también
Categorías
Más información sobre Dates and Time en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!