my code doesnt run because the sizes are different

3 Ansichten (letzte 30 Tage)
Vasco
Vasco am 9 Sep. 2023
Beantwortet: Vasco am 9 Sep. 2023
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 Kommentare
Torsten
Torsten am 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.
Vasco
Vasco am 9 Sep. 2023
the error is this one
Unable to perform assignment because the size of the left side is 360-by-90 and the size of the right side is 360-by-180.
Ill try to adapt the code but im almost imploding i have 1 day to send a report and im handling the methodology as im doing the code and im handling mistakes and improving the writing il try to do it fast sorry.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Vasco
Vasco am 9 Sep. 2023
This is the code that I managed to get
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/2, ntime);
for i = 1:length(lista_ficheiros)
filename = fullfile(input, lista_ficheiros(i).name);
nc = netcdf.open(filename);
lwe = ncread(filename, 'lwe_thickness');
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);
lat_idx = lat >= 0 & lat <= 90;
lat = lat(lat_idx);
lwe = lwe(:, 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, :, :));
contourf(lon, lat, data', clevv, 'LineStyle', 'none', 'Fill', 'on'); % Transpose data
xlabel('Longitude');
ylabel('Latitude');
title('Mapa da média de lwe');
clim([min(clevv), max(clevv)]);
colormap(hot(length(clevv) - 1));
colorbar('eastoutside');
c = colorbar;
c.Label.String = 'Valor em metros (m)';
ylim([0 90]);
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);
clevvv = 0:0.01:0.2;
contourf(lon, lat, pixel_std', clevvv, 'LineStyle', 'none', 'Fill', 'on'); % Transpose final2
xlabel('Longitude');
ylabel('Latitude');
clim([min(clevvv), max(clevvv)]);
colormap(hot(length(clevvv) - 1));
colorbar('eastoutside');
c = colorbar;
c.Label.String = 'Valor em metros (m)';
title('Standard deviation');
hAx = gca;
hAx.YDir = 'normal';
ylim([0 90]);
hold on
plot(C.coastlon, C.coastlat, 'k')
figura_4 = fullfile(output, strcat(nome_std, '.png'));
saveas(gcf, figura_4);

Weitere Antworten (1)

Walter Roberson
Walter Roberson am 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 Kommentare
Vasco
Vasco am 9 Sep. 2023
the 6 months of data i mean the year 2017 only has 5 files which only equate to 6 months of data
Vasco
Vasco am 9 Sep. 2023
UPDATE
I managed to do it

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Write Data to Channel finden Sie in Help Center und File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by