weighted variable on a satellite

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);

Index exceeds the number of array elements. Index must not exceed 0.

Answer 1

[Walter Roberson]

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.

[Vasco]

that part doesnt run into any problems, i dont have any problems on running I understand that it may not be optimized but its the best i can do with a limited time

What I cant do is the cosd(lat) * by the lwe_data to take into account the corrections of the pixels.

please guide me on it

[Walter Roberson]

You have an array lwe which is length(lon) by length(lat) . You want to multiply each element in the array by the cosd() of its corresponding latitude, cosd(lat) . You are not wanting to do algebraic multiplication, you want element-by-element multiplication, so you will be wanting to use .* cosd(LAT) instead of * cosd(LAT) where LAT is element-by-element the latitude that applies for that location. Which... you already calculated and placed into the y variable.

lwe .* cosd(y)

[Vasco]

I came up with this code but the values that are appearing are to big

cos_lat = cosd(lat);
lwe_peso = lwe .*cos_lat';

nvall = nnz(~isnan(lwe_peso));
ss = sum(lwe_peso(~isnan(lwe_peso)));
media_s = sum(ss(~isnan(ss))) / nvall ; 
total = sum(cos_lat(~isnan(cos_lat)));
globalm = ss / total; 

med_global(i) = globalm;

[Vasco]

in a temporal series my values where about 0.09 and now theu show as 9 and i dont know that it means

does the unit changed all the variables appear in metres it could only be explained if it turned into centimeters and it couldnt have happened

[Walter Roberson]

that part doesnt run into any problems

You are mistaken on that point. Running the right calculation on the wrong data files is no virtue.

Example filename: GRD-3_2002094-2002120_GRAC_JPLEM_BA01_0600_LND_v04.nc

Y = input('year to examine');
SD = input('start day in year');
ED = input('end day in year');

input_folder = 'C:\Users\vasco\OneDrive\Ambiente de Trabalho\ESTAGIO_Vasco\JPL TELLUS GRACE Level-3\';
output_folder = 'C:\Users\vasco\OneDrive\Ambiente de Trabalho\ESTAGIO_Vasco\JPL TELLUS GRACE Level-3 Output\';
lista_ficheiros = dir(fullfile(input_folder, '*.nc'));

filenames = {lasta_ficherios.name};
name_pieces = regexp(filenames, '_', 'split');
date_ranges = cellfun(@(NP) NP(2), name_pieces);
BYs = cellfun(@(NP) str2double(NP(1:4)), date_ranges);
EYs = cellfun(@(NP) str2double(NP(9:12)), date_ranges);
SDs = cellfun(@(NP) str2double(NP(5:6)), date_ranges);
EDs = cellfun(@(NP) str2double(NP(13:14)), date_ranges);
BD = datetime(BYs, 1, 0) + days(SDs);
ED = datetime(EYs, 1, 0) + days(EDs);

Now BD and ED are datetime arrays. BD is the starting date that each file pertains to, and ED is the ending date that each file pertains to. You can use isbetween to determine whether you want to pay attention to the file.

Because each file is for a range of dates, you have the problem that you might be wanting to focus on a range of days that starts after the beginning of a file, or that ends before the end of a file. Reading only part of a .nc file can take more work, depending on how the file is structured.

Or... what you could be doing is using ismember() on a file name to determine whether the file name is one you wanted to exclude...

Or... what you could do is use uigetfile() with 'MultiSelect', 'on', and have the user select which nc files to process.

[Vasco]

Ill try to fix it but in the mean time can you help me do the code that regards this:

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.

I this is gonna be the second night i havent slept i need something to justify this mindless work hours to myself

[Vasco]

this is the only thing i need to finish my college degree.

[Walter Roberson]

ss = sum(lwe_peso(~isnan(lwe_peso)));
media_s = sum(ss(~isnan(ss))) / nvall ; 
total = sum(cos_lat(~isnan(cos_lat)));

isnan(lwe_peso) is going to be a 2d logical array. When you index an array of any dimension at a 2d logical array as a single subscript, the result is a column vector, never a 2D array. So lwe_peso(~isnan(lwe_peso)) would be a column ector and sum() of that would be a scalar, so ss is going to be a scalar.

You selected only non-nan elements of lwe_peso so the only way that the sum() of those non-nan elements could possibly be nan, would be if there were infinities of opposite sign. inf + (-inf) --> nan . If you have infinities in the data you should probably be treating them the same as a nan. In particular you should consider using lwe_peso(isfinite(lwe_peso)) instead of lwe_peso(~isnan(lwe_peso)) .

Since the result of the lwe_peso indexing is always going to be a column vector, then instead of bothering dong the non-nan selecting of the data yourself, you should consider using

ss = sum(lwe_peso, 'all', 'omitnan');

For the sample .nc file you posted, the entries are all finite so the (scalar) sum will be finite, so there is no point taking sum(ss(~isnan(ss)) -- isnan will never be true (if lwe contained only finite values) so you would effectively being taking sum(ss) where ss is already a scalar.

If you are not doing any further calculations with ss you should consider replacing those series of lines with

media_s = mean(lwe_peso, 'all', 'omitnan');

[Vasco]

i checked the lwe_peso in the workspace and the calculations for multiplying by cos_lat are wrong i truly dont know what i can possibly do to overcome this to be honest with you

thank you for your time but im in a stand still with this

[Vasco]

update

all this time the problem that i was having was because of the fact that the lat value is alwyas .5 and not integer, i am very sorry for wasting your time

my deepest apologies

now i am presented with a new challenge