Collect data from multiple sets of lat and long instead of just one set

I have this code below in matlab, it collects data from AOD in a set of lat and long, however, I want to collect in multiple locations of lat and lon. Do you have any idea? I have this code below in matlab, it collects data from AOD in a set of lat and long, however, I want to collect in several places of lat and lon. Do you have any idea? The script I used is below.

The data looks like this: MYD04_3K.A2010001.1540.061.2018053201411.hdf

%%
% MYD04 - DADOS PARTE 1

% DATAS -> MYD04_3K.A2010001.1540.061.2018053201411.hdf

%% PROGRAMA_LER_HDF_SÉRIE_TEMPORAL_AOD_MODIS

%  OBETIVO: ler arquivos HDF e extrair a série temporal da variável desejada

%  DESCRIÇÃO: abre um arquivo HDF de um diretório específico,identifica dentro do HDF as variáveis que se deseja trabalhar 
%  (ex.: Latitude, Longitude, optical_depth_land_and_ocean), escreve esses arquivos de forma(double), identifica nesses arquivos 
%  os valores de erro/desprezíveis (ex.: -9999), substitui por "NaN", gera um arquivo corrigido, multiplica os dados dessa matriz 
%  pelo fator de conversão especificado no arquivo HDF (ex.:0.01), gera um novo arquivo(com os valores físicos reais) e em seguida plota os gráficos

%% REGISTRA_DIRETÓRIOS

clc
clear all
% close all
%%
% DirDSc=pwd;

DirMYD = '/media/augusto/SATELITES/MYD043K'
% DirDMY='D:\MODIS-ATM\MODIS043K\MYD043K';
% DirDMI='D:\MODIS-ATM\MCD19A2\MCD19A2d';

% LOCALIZA_ÁREAS_ESTUDO
% Encontra os pontos de lat e lon 

%Localizações
%Bel: -1.473137 -48.485657
%Stm: -2.452059 -54.753617
%Alt: -3.217477 -52.230464
%Mab: -5.394255 -49.135187
%ParaC: -4.27 -52.60

DirShp='/media/augusto/AUGUSTOG/IBGEShp/bcim_2016_shapefiles_21-11-2018';cd(DirShp);
shpc = shaperead('loc_capital_p.shp'); %capitais AS

cgeo={shpc.X;shpc.Y;shpc.nome}';
disp(cgeo);

cgeo=input('Informe a coordenada (graus) do local de estudo /Ex: [-02.133333 -58.983609]:');
pixe_l=input('Informe o pixel, cujo centro é local de estudo /Ex: [15]:');

lat_n=cgeo(1,1)+((pixe_l(1,1)./2)./111.1);
lat_s=cgeo(1,1)-((pixe_l(1,1)./2)./111.1);

lon_e=cgeo(1,2)+((pixe_l(1,1)./2)./111.1);
lon_w=cgeo(1,2)-((pixe_l(1,1)./2)./111.1);

%%
% LER_HDF

% Ler o arquivo HDF em um diretório específico. 
% Ler as variáveis desejadas do grupo de variáveis disponibilizadas no arquivo HDF que foi aberto.

% 1) Latitude
% 2) Longitude
% 3) Optical_Deph_Land_and_Ocean
% 4) Solar_Zenith

t1=datestr(now);
disp(t1)

cd(DirMYD); %Terra

m=dir('*hdf*'); %diretorio das pastas

MYD04L2=[];
MYD04L2e=[];

% MYD04L2=[];
% MYD04L2e=[];
%%
for j=1:length(m)
    jj=num2str(j);
    
    OpenFile{1,1}=num2str(j);
    OpenFile{1,2}=m(j).name;
    disp(OpenFile) %mostra na tela posição processamento
        
    ano(j)= str2double(m(j).name(11:14));     %#ok<SAGROW>
    dia(j)= str2double(m(j).name(15:17));  %#ok<SAGROW>

    %   MOD_DAY=[];for i=3:length(dir);MOD_DAY=[MOD_DAY;[str2num(f(i).name(19:22))]];end;    
    %--------------------------------------------------------------------------    

    % Especificando a precisão dos números /nomeando variáveis HDF na workspace

        cmd=['lat = double(hdfread(''' m(j).name ''', ''/mod04/Geolocation Fields/Latitude''));'];eval(cmd);
        cmd=['lon = double(hdfread(''' m(j).name ''', ''/mod04/Geolocation Fields/Longitude''));'];eval(cmd);
        cmd=['asz = double(hdfread(''' m(j).name ''', ''/mod04/Data Fields/Solar_Zenith''));'];eval(cmd);
        cmd=['aod = double(hdfread(''' m(j).name ''', ''/mod04/Data Fields/Optical_Depth_Land_And_Ocean''));'];eval(cmd);

    %--------------------------------------------------------------------------    

    % Substituindo erro p/ NaN    
    % Encontrando os índices (posição nas matrizes que possuem valores negativos e/ou menores que -9999)

    id_lat=find(lat == -9999);
    id_lon=find(lon == -9999);
    lat(id_lat)=NaN;  %#ok<SAGROW>
    lon(id_lon)=NaN; %#ok<SAGROW>

    id_aod=find(aod<=0);
    aod(id_aod)=NaN;   %#ok<SAGROW>
    aod=0.001.*aod;

    id_asz=find(asz<=-9999);
    asz(id_asz)=NaN; %#ok<SAGROW>
    asz=0.01.*asz;

    %--------------------------------------------------------------------------

    %Selecionando_pixel e a profundidade óptica média (representativa do pixel) 
    a_pixe=find(lat>=lat_s & lat<=lat_n & lon>=lon_w & lon<=lon_e);
    
    cmd=['Dnum=[datenum(''31-Dez-' num2str(ano(j)-1) ' 00:00:00'') + (str2double(m(j).name(15:17)) + str2double(m(j).name(19:20))./24 + str2double(m(j).name(21:22))./60./24)];'];eval(cmd); % Date(n)
    Dvec=datevec(Dnum); %DateVec
    JDay=(str2double(m(j).name(15:17)) + str2double(m(j).name(19:20))./24 + str2double(m(j).name(21:22))./60./24); %JulianDate

    MYD04L2=[MYD04L2;[Dnum,Dvec,JDay,nanmean(aod(a_pixe))]];  % p/TERRA usar MOD04L2
    MYD04L2e=[MYD04L2e;[NaN(size(a_pixe)),lat(a_pixe),lon(a_pixe),(aod(a_pixe))]];  % p/TERRA usar MOD04L2
    MYD04L2e(isnan(MYD04L2e(:,1)),1)=Dnum; 
    
%     MYD04L2=[MYD04L2;[Dnum,Dvec,JDay,nanmean(aod(a_pixe))]];  % p/AQUA usar MYD04L2
%     MYD04L2e=[MYD04L2e;[NaN(size(a_pixe)),lat(a_pixe),lon(a_pixe),(aod(a_pixe))]];  % p/AQUA usar MYD04L2
%     MYD04L2e(isnan(MYD04L2e(:,1)),1)=Dnum;    
    
end
%%
cd ..
clear('id_aod','id_asz','ans','cmd','id_lat','id_lon');

t2=datestr(now);

%% MATRIZ FINAL DADOS MODIS (MOD04L23K)

MYDIS04L2=[MYD04L2]% MODIS04L2=[MOD04L2;MYD04L2]; % empilhando as matrizes
[ord,ind]=sort(MYDIS04L2(:,1)); % ordenando as medidas
MYDIS04L2=MYDIS04L2(ind,:);

%plot(MYDIS04L2(:,1),MYDIS04L2(:,9),'.')

[Augusto Gabriel da Costa Pereira]

I tried to look for a solution through chatgpt and the following was suggested:

Swap these lines below:

cgeo=input('Informe a coordenada (graus) do local de estudo /Ex: [-02.133333 -58.983609]:');
pixe_l=input('Informe o pixel, cujo centro é local de estudo /Ex: [15]:');

lat_n=cgeo(1,1)+((pixe_l(1,1)./2)./111.1);
lat_s=cgeo(1,1)-((pixe_l(1,1)./2)./111.1);

lon_e=cgeo(1,2)+((pixe_l(1,1)./2)./111.1);
lon_w=cgeo(1,2)-((pixe_l(1,1)./2)./111.1);

And insert these:

coords = [-1.473137 -48.485657;
          -2.452059 -54.753617;
          -3.217477 -52.230464;
          -5.394255 -49.135187;
          -4.27 -52.60;
          -1.2345 -45.6789]; % Adicione aqui as 6 latitudes e longitudes específicas que você deseja coletar

pixe_l=input('Informe o pixel, cujo centro é local de estudo /Ex: [15]:');

for i = 1:size(coords, 1)
    lat_n(i) = coords(i, 1) + ((pixe_l(1,1)./2)./111.1);
    lat_s(i) = coords(i, 1) - ((pixe_l(1,1)./2)./111.1);
    lon_e(i) = coords(i, 2) + ((pixe_l(1,1)./2)./111.1);
    lon_w(i) = coords(i, 2) - ((pixe_l(1,1)./2)./111.1);
end

But I was not successful

[dpb]

As per usual, it's essentially impossible to do anything data specific such as your request without an actual file to work on. Attach a representative file so somebody here can access it.

Also give a specific example of what it is you would like to return from the file, generalities don't really define a a problem sufficiently.

It would appear the suggested code above in the loop is defining a specific range of coordinates, but there's no indication in your code that you then did anything with those inside that loop or elsewhere with a loop over than list of coordinates to actually use them to return any data???

Unless the data file is too large to hold in memory, the easiest approach might be to read the whole file and then just use logical addressing to save the sections of interest.

[Stephen23]

Get rid of all anti-pattern EVAL calls. For example, replace this:

cmd=['lat = double(hdfread(''' m(j).name ''', ''/mod04/Geolocation Fields/Latitude''));'];

eval(cmd);

with simpler, neater, more robust, much more efficient:

lat = double(hdfread(m(j).name, '/mod04/Geolocation Fields/Latitude'));

and then replace all of those deprecated DATENUM/DATEVEC/DATESTR with DATETIME/DURATION. You can plot them directly too.

Code like this is hard to write and even harder to maintain:

cmd=['Dnum=[datenum(''31-Dez-' num2str(ano(j)-1) ' 00:00:00'') + (str2double(m(j).name(15:17)) + str2double(m(j).name(19:20))./24 + str2double(m(j).name(21:22))./60./24)];'];

eval(cmd);

The intent is rather obfuscated, so that I will need some time to reverse-engineer it and figure out what you are actually doing. ... EDIT: okay, it appears that these are date/time units in the filenames:

m(j).name(11:14) % year

m(j).name(15:17) % doy

m(j).name(19:20) % hours

m(j).name(21:22) % minutes

Lets try it now with your example filename:

j = 1;
m.name = 'MYD04_3K.A2010001.1540.061.2018053201411.hdf';
% This one line is all you need:
dt = datetime(m(j).name(11:22), 'InputFormat','uuuuDDD.HHmm'); % simple!
% Just to compare:
dt.Format = 'u MM dd HH mm ss'

dt = datetime
   2010 01 01 15 40 00

vs. that complex, obfuscated, fragile (it only works in your locale!), hard-to-debug, inefficient code:

ano(j)= str2double(m(j).name(11:14))

ano = 2010

cmd=['Dnum=[datenum(''31-Dec-' num2str(ano(j)-1) ' 00:00:00'') + (str2double(m(j).name(15:17)) + str2double(m(j).name(19:20))./24 + str2double(m(j).name(21:22))./60./24)];'];
eval(cmd); %             ^^^ ugh, locale dependent!
% Just to compare:
dv = datevec(Dnum)

dv = 1×6

2010 1 1 15 40 0

So in half the code I wrote a more robust date conversion using the recommended DATETIME. You should too.

A few other improvements:

Replace CD (slow and hard to debug) with absolute/relative filenames using FULLFILE.

Replace CLEAR ALL with CLEARVARS (why do you need to clear cached functions?)

I get the feeling you are doing some kind of interpolation or something similar, but it is rather obfsucated. Refactoring is probably a good idea.

[dpb]

'Dnum=[datenum(''31-Dez-' num2str(ano(j)-1) ' 00:00:00'') + ...' ...
               '(str2double(m(j).name(15:17)) + ...' ...
               ' str2double(m(j).name(19:20))./24 + ...' ...
               ' str2double(m(j).name(21:22))./60./24)];'];

Well, breaking it into pieces...

datenum(''31-Dez-' num2str(ano(j)-1) ' 00:00:00'') 

is the day for the year one less than that in the array, so this could be

Dnum=datenum(ano(j)-1),12,31);

replacing the need to parse the string with the equivalent month, day values...

The rest must be adding time to the day since a datenum time is the fractional part of the 24-hr day...let's see what the next piece is -- oh! the year-1 is that the filename is day of year, hour, minute it appears so crudely the above would translate to

da=str2double(m(j).name(15:17));
hr=str2double(m(j).name(19:20))/24;
mn=str2double(m(j).name(21:22))/60/24;
Dnum=datenum(ano(j)-1),12,31+da,hr,mn);

I'm guessing this would also be

da=str2double(m(j).name(15:17));
hr=str2double(m(j).name(19:20))/24;
mn=str2double(m(j).name(21:22))/60/24;
Dnum=datenum(ano(j)),1,1+da-1,hr,mn);

which more simply could surely be written as if the file name happened to include the above day and time as

>> m.name(15:22)='003_0230';        % make up a partial file name with the date portion
>> yr=string(ano(1));               % and extract the year from the array, turn to string
DT=datetime(yr+m.name(15:22),'InputFormat','uuuuDDD_HHmm') % combine year, day-time and convert
DT = 
  datetime
   03-Jan-2022 02:30:00
>> 

See <datetime> and follow links to input formats for the formatting string conventions; they are NOT the same as for datenum.

[Augusto Gabriel da Costa Pereira]

@Stephen23

I see no need to make such adjustments now. I'm relatively new to the programming world, so there's no need for you to write words in a mocking and derogatory tone towards my script.

For a beginner I'm pretty good, when you were starting out you certainly weren't the super smart one.

So be careful with your words, the script is working perfectly, I hope you're not a teacher, if you are, your students must be sadder precisely because of your sarcastic and demotivating comments!!!

I didn't ask for a solution to what you suggested, I just asked for a solution to adapt the script to perform AOD collection for multiple lat and long sets and not just one as it is at the moment. If you have a real idea to add, submit it, otherwise be silent.

Indeed, you are a poet when you are silent.

[Augusto Gabriel da Costa Pereira]

@Mathieu NOE, @Image Analyst and @Chad Greene, do you have a solution?

The data used is attached to Google Drive, see the link below: https://drive.google.com/drive/folders/18GgPs-iN3brZHQ3pZNK8GyLt4SOmemY8?usp=sharing

[Augusto Gabriel da Costa Pereira]

@dpb this link is here: https://drive.google.com/drive/folders/18GgPs-iN3brZHQ3pZNK8GyLt4SOmemY8?usp=sharing

[Mathieu NOE]

hello @Augusto Gabriel da Costa Pereira

sorry I cannot run the code (I don't have the Mapping Toolbox)

[dpb]

Nor I...