How to plot satellite data?

I am currently working with IDL to plot spectral images of trace gases concentration from satellite observation (such as TROPOMI on S5P and IASI on Metop satellites) but I want to switch to MATLAB. In this process I am having several trouble:

1. I can not find a function or a tool to read native format images (.nat extension).
2. I read files in NetCDF or HDF5 formats but I can not find a good way to plot this data.

I attacched an example of what I am trying to obtain. If some one has an sample code I would be grateful if he shared it.

[Clare Davis]

Hello, I am having the same trouble - did you get a good answer in the end? Thanks

[Meg Noah]

I do not have a .nat extension file. Can you post one here?

Answer 1

[Meg Noah]

Here's a solution

clc
close all
clear all

filename = 'W_XX-EUMETSAT-Darmstadt,HYPERSPECT+SOUNDING,METOPA+IASI_C_EUMP_20181108012658_62544_eps_o_l2.nc';
info = ncinfo(filename);

% to get a listing on console of variable names
varnames =  {info.Variables.Name}';
idx = find(contains(varnames,'integrated_co') & ~contains(varnames,'integrated_co2'));

% learning about the data
Attributes = info.Variables(idx).Attributes;
for i = 1:length(Attributes)
    if (isnumeric(Attributes(i).Value))
    fprintf(1,'%s = %f\n',Attributes(i).Name,Attributes(i).Value);
    else
    fprintf(1,'%s = %s\n',Attributes(i).Name,Attributes(i).Value);
    end
end

dataCO = ncread(filename,'integrated_co');
% they seem to have made the missing data value a really big value
minval = min(dataCO(:));
dataCO(dataCO == max(dataCO(:))) = -300;
maxval = max(dataCO(:));
fprintf(1,'minval = %f\n',minval);
fprintf(1,'maxval = %f\n',maxval);
minval = 2e-4;
maxval = 67e-4;
m = 63/(maxval-minval);
b = -m*minval;
scaleddataCO = uint8(m*dataCO+b);
scaleddataCO(dataCO<0) = nan; % missing data
colorvals = colormap(jet);
lat = ncread(filename,'lat');
lon = ncread(filename,'lon');
[ny,nx] = size(dataCO);
figure()
hold on
for i = 1:64
    idx = find(scaleddataCO == i);
    plot(lon(idx),lat(idx),'s', 'markerfacecolor',colorvals(i,:),'markeredgecolor','none');
end
xlim([-180 180])
ylim([-90 90]);
axis equal
axis tight
set(gca,'ytick',[-90 -60 -30 0 30 60 90]);
set(gca,'xtick',[-180 -150 -120 -90 -60 -30 0 30 60 90 120 150 180]);
grid on
box on
title('Integrated CO 10^{-4} kg/m^2');
colormap(jet)
hb = colorbar;
hb.Ticks = [0 0.2 0.4 0.6 0.8 1];
vals = linspace(2,67,6);
for j = 1:length(vals)
hb.TickLabels{j} = num2str(vals(j));
end

Here's how to plot it with map functions

hfig = figure('color','white','position',[32 200 1200 650]);
worldmap('World')
load coastlines
plotm(coastlat,coastlon,'k')
hfig.Name = 'Integrated CO';
hold on;
for i = 1:64
    idx = find(scaleddataCO == i);
    if (~isempty(idx))
    plotm(lat(idx),lon(idx),'s', 'markerfacecolor',colorvals(i,:),'markeredgecolor','none');
    end
end
title('Integrated CO 10^{-4} kg/m^2');
colormap(jet)
hb = colorbar('southoutside');
hb.Ticks = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1];
vals = linspace(2,67,11);
for j = 1:length(vals)
hb.TickLabels{j} = num2str(vals(j));
end

[Visweshwaran R]

@Meg Noah, the code is well written and precise. I too faced the same problem in plotting the netCDF file for soil moisture data (ASCAT). While I try to adjust this code for my datasets, I faced a few problems.

1) After printing the min and max value. you have assigned the same with 2e-4 and 67e-4.What is the need for it?

2)While calculating 'm', you used 63, why?

3)Similarly, you used for loop till 64, why that number particularly?

I know these questions may be basic but I am struck and 5 min of your time will be really helpful to me.

Thank you in advance.

[Hugo]

What would it look like without the missing data part?

I'm having trouble adapting the code to a random satellite ncDEF file cause it doesnt plot all my data

Answer 2

[KSSV]

Reading netCDF files is very easy in MATLAB. You need to just know ncdisp, ncread ..once you have read the data use pcolor to get the data you wanted to plot.

[A.G. Bruno]

My problem is that my data are 3 matrix of the same dimension, 2 for the coordinates of the satellite track and one for data observed.

[KSSV]

Let it be......attach your file....

[A.G. Bruno]

clear all
close all

yyyy = '2018';
mm = '11';
dd = '08';

folder = '...'
file1 = '...'

filename = dir(strcat(folder, file1))

data = strcat(folder, filename.name)
    
[finfo outstrct]=read_nc_file_struct(data)
    
lon = outstrct.lon;
lat = outstrct.lat;
int_co = outstrct.integrated_co;
int_co = int_co';

This code is for a single file, to read .nc files I am currently try to use the routine read_nc_file_struct suggested by Sentinel staff. The output structures have the following dimensions.

outstrct = 

  struct with fields:

                          filename: '/scratch/atmcomp/shared/iasi/L2_metopa_nc/2018/11/08/W_XX-EUMETSAT-Darmstadt,HYPERSPECT+SOUNDING,METOPA+IASI_C_EUMP_20181108232058_62557_eps_o_l2.nc'
                   cloud_formation: [3×1 int16]
              pressure_levels_temp: [101×1 double]
          pressure_levels_humidity: [101×1 double]
             pressure_levels_ozone: [101×1 double]
    surface_emissivity_wavelengths: [12×1 double]
                 record_start_time: [765×1 double]
                  record_stop_time: [765×1 double]
                  degraded_ins_MDR: [765×1 double]
                 degraded_proc_MDR: [765×1 double]
                               lat: [120×765 double]
                               lon: [120×765 double]
                      solar_zenith: [120×765 double]
                  satellite_zenith: [120×765 double]
                     solar_azimuth: [120×765 double]
                 satellite_azimuth: [120×765 double]
        fg_atmospheric_temperature: [101×120×765 double]
        fg_atmospheric_water_vapor: [101×120×765 double]
              fg_atmospheric_ozone: [101×120×765 double]
            fg_surface_temperature: [120×765 double]
           atmospheric_temperature: [101×120×765 double]
           atmospheric_water_vapor: [101×120×765 double]
                 atmospheric_ozone: [101×120×765 double]
               surface_temperature: [120×765 double]
            integrated_water_vapor: [120×765 double]
                  integrated_ozone: [120×765 double]
                    integrated_n2o: [120×765 double]
                     integrated_co: [120×765 double]
                    integrated_ch4: [120×765 double]
                    integrated_co2: [120×765 double]
                surface_emissivity: [12×120×765 double]
           number_cloud_formations: [120×765 double]
            fractional_cloud_cover: [3×120×765 double]
             cloud_top_temperature: [3×120×765 double]
                cloud_top_pressure: [3×120×765 double]
                       cloud_phase: [3×120×765 double]
                  surface_pressure: [120×765 double]
                   instrument_mode: [765×1 double]
               spacecraft_altitude: [765×1 double]
                      flag_amsubad: [120×765 double]
                     flag_avhrrbad: [120×765 double]
                       flag_cdlfrm: [120×765 double]
                       flag_cldnes: [120×765 double]
                       flag_cdltst: [120×765 uint16]
                       flag_daynit: [120×765 double]
                      flag_dustcld: [120×765 double]
                      flag_fgcheck: [120×765 uint16]
                      flag_iasibad: [120×765 double]
                       flag_initia: [120×765 double]
                       flag_itconv: [120×765 double]
                      flag_landsea: [120×765 double]
                       flag_mhsbad: [120×765 double]
                        flag_numit: [120×765 double]
                       flag_nwpbad: [120×765 double]
                    flag_physcheck: [120×765 double]
                     flag_retcheck: [120×765 uint16]
                       flag_satman: [120×765 double]
                      flag_sunglnt: [120×765 double]
                       flag_thicir: [120×765 double]
                       nerr_values: [765×1 double]
                  error_data_index: [120×765 double]
                 temperature_error: [406×116×765 uint32]
                water_vapour_error: [171×116×765 uint32]
                       ozone_error: [55×116×765 uint32]
                         surface_z: [120×765 double]
                          co_qflag: [120×765 double]
                           co_bdiv: [120×765 uint32]
                           co_npca: [120×765 double]
                     co_nfitlayers: [120×765 double]
                     co_nbr_values: [765×1 double]
                         co_cp_air: [19×120×765 double]
                        co_cp_co_a: [19×120×765 double]
                           co_x_co: [19×120×765 double]
                  co_h_eigenvalues: [10×120×765 double]
                 co_h_eigenvectors: [190×120×765 double]

[KSSV]

pcolor(lon,lat,int_co) ;
shading interp 
colorbar

This doesnt give you a plot?

[A.G. Bruno]

No, the problem is the dimension of the arrays.

Error using pcolor (line 61)
Matrix dimensions must agree.

Error in prova20122018_v2 (line 53)
pcolor(lon,lat,int_co) ;

[KSSV]

What are the dimensions of each variable?

Try:

pcolor(lon,lat,int_co') ;

I have asked you to attach the data/ file so that you can get more help straight away.

[A.G. Bruno]

It still doesn't work.

Sorry, I did not understand that you mean the data file. I've attached a dropbox link because the file is too large to upload it.

https://www.dropbox.com/s/2dkep0bzimms66m/W_XX-EUMETSAT-Darmstadt%2CHYPERSPECT%2BSOUNDING%2CMETOPA%2BIASI_C_EUMP_20181108012658_62544_eps_o_l2.7z?dl=0

Answer 3

[ming lei]

I also encountered a similar problem. Is it finally solved?