my code doesnt run because the sizes are different
3 views (last 30 days)
Show older comments
Im running a code where I want to study the northern hemishpere but the code doesnt run and i cant fix please help me ill put the code and one of the files so you can see it.
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');
idx = lon>180;
lon(idx) = lon(idx) - 360;
avanco = 1-find(idx,1);
lon = circshift(lon,avanco);
lwe = circshift(lwe,avanco,1);
lon_idx = lon>=-180 & lon<=180;
lat_idx = lat>0 & lat<90;
lon = lon(lon_idx);
lat = lat(lat_idx);
lwe = lwe(lon_idx,lat_idx);
lwe_data(:, :, i)=lwe;
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_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);%%
2 Comments
Torsten
on 9 Sep 2023
You want an answer ? Then make it easy for us to handle your problem.
Please adapt your code and the input file so that it can be run without further changes from our side to reproduce the error you encounter.
Accepted Answer
More Answers (1)
Walter Roberson
on 9 Sep 2023
This code will fail when it gets to the plotting section, as the plotting section assumes that at least 158 files have been read in. You should find some other way of selecting a subset of the data.
format long g;
inputfolder = '.';
outputfolder = '.';
lista_ficheiros = dir(fullfile(inputfolder, '*.nc'));
lon_1 = double(ncread(fullfile(inputfolder, lista_ficheiros(1).name), 'lon'));
lat_1 = double(ncread(fullfile(inputfolder, lista_ficheiros(1).name), 'lat'));
nlon = length(lon_1);
nlat = length(lat_1);
ntime = length(lista_ficheiros);
for i = 1:length(lista_ficheiros)
filename = fullfile(inputfolder, 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');
idx = lon>180;
lon(idx) = lon(idx) - 360;
avanco = 1-find(idx,1);
lon = circshift(lon,avanco);
lwe = circshift(lwe,avanco,1);
lon_idx = lon>=-180 & lon<=180;
lat_idx = lat>0 & lat<90;
lon = lon(lon_idx);
lat = lat(lat_idx);
lwe = lwe(lon_idx,lat_idx);
if i == 1
lwe_data = zeros([size(lwe), length(lista_ficheiros)]);
end
lwe_data(:, :, i)=lwe;
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(outputfolder, strcat(nome_ficheiro, '.png'));
saveas(gcf, figura);
estatisticas = fullfile(outputfolder, 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(outputfolder, 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(outputfolder, '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_outputs = fullfile(outputfolder, '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(outputfolder, 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(outputfolder, strcat(nome_std, '.png'));
saveas(gcf, figura_4);%%
6 Comments
See Also
Categories
Find more on Write Data to Channel in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!