How to extract data with latitude and longitude?

17 views (last 30 days)
Hello, I used this code to grid swath OMI data. 1 swath file have data with size: 60x1643 (we have 15 files swath per day) and after run this code the globle size: 60x24856.matrix1.JPG
and I got the picture of globle:
OMI-Aura_L2-OMNO2_2014m1231t2351-o55654_v003-2018m0617t035839.he5.m.png
This is the code:
thepath='D:\Paper_NOx\newtest';
% Open the HDF5 File.
FILE_NAME = ...
'OMI-Aura_L2-OMNO2_2014m1231t2351-o55654_v003-2018m0617t035839.he5';
file_id = H5F.open (FILE_NAME, 'H5F_ACC_RDONLY', 'H5P_DEFAULT');
% Open the dataset.
DATAFIELD_NAME = '/HDFEOS/SWATHS/ColumnAmountNO2/Data Fields/ColumnAmountNO2';
Lat_NAME='HDFEOS/SWATHS/ColumnAmountNO2/Geolocation Fields/Latitude';
Lon_NAME='HDFEOS/SWATHS/ColumnAmountNO2/Geolocation Fields/Longitude';
data_id = H5D.open (file_id, DATAFIELD_NAME);
% Read the units.
ATTRIBUTE = 'Units';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
units = H5A.read(attr_id, 'H5ML_DEFAULT');
% Read the offset.
ATTRIBUTE = 'Offset';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
offset = H5A.read(attr_id, 'H5ML_DEFAULT');
% Read the scale.
ATTRIBUTE = 'ScaleFactor';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
scale = H5A.read(attr_id, 'H5ML_DEFAULT');
% Read the fill value.
ATTRIBUTE = '_FillValue';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
fillvalue=H5A.read (attr_id, 'H5T_NATIVE_DOUBLE');
% Read the missing value.
ATTRIBUTE = 'MissingValue';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
missingvalue=H5A.read (attr_id, 'H5T_NATIVE_DOUBLE');
% Read title attribute.
ATTRIBUTE = 'Title';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
long_name=H5A.read (attr_id, 'H5ML_DEFAULT');
% Close and release resources.
H5A.close (attr_id)
H5D.close (data_id);
H5F.close (file_id);
f = figure('Name','OMI_NO2_merged','visible','off');
% Create the plot.
axesm('MapProjection','eqdcylin',...
'Frame','on','Grid','on', ...
'MeridianLabel','on','ParallelLabel','on','MLabelParallel','south')
coast = load('coast.mat');
D = dir(fullfile(thepath, '*.he5'));
for k =1:numel(D)
FILE_NAME1 = fullfile(thepath, D(k).name);
file_id1 = H5F.open (FILE_NAME1, 'H5F_ACC_RDONLY', 'H5P_DEFAULT');
data_id1 = H5D.open (file_id1, DATAFIELD_NAME);
lat_id=H5D.open(file_id1, Lat_NAME);
lon_id=H5D.open(file_id1, Lon_NAME);
% Read the dataset.
data=H5D.read (data_id1,'H5T_NATIVE_DOUBLE', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT');
lat=H5D.read(lat_id,'H5T_NATIVE_DOUBLE', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT');
lon=H5D.read(lon_id,'H5T_NATIVE_DOUBLE', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT');
% Convert the data to double type for plot.
data=double(data);
lon=double(lon);
lat=double(lat);
% Replace the filled value with NaN.
data(data==fillvalue) = NaN;
% Multiply scale and adding offset, the equation is scale *(data-offset).
data = scale*(data-offset);
% surfm is faster than contourfm.
surfm(lat, lon, data);
if k == 1
lat_m = [lat];
lon_m = [lon];
data_m = [data];
else
lat_m = [lat_m, lat];
lon_m = [lon_m, lon];
data_m = [data_m, data];
end
% Detach from the Swath object.
H5D.close (data_id1);
H5F.close(file_id1);
end
colormap('Jet');
h=colorbar();
plotm(coast.lat,coast.long,'k')
% Draw unit.
set(get(h, 'title'), 'string', 'None', ...
'FontSize', 12, 'FontWeight','bold', ...
'Interpreter', 'none');
% Put title.
tstring = {'OMI_Aura_L2 Merged';DATAFIELD_NAME};
title(tstring, 'Interpreter', 'none', 'FontSize', 16, ...
'FontWeight','bold');
% The following fixed-size screen size will look better in JPEG if
% your screen is too large.
scrsz = [1 1 800 600];
set(f, 'position', scrsz, 'PaperPositionMode', 'auto');
saveas(f, [FILE_NAME '.m.png']);
Now I want to extract data of USA region but I don't know how to extract data with lat and lon. Please give me solution to solves it. Thank you so much!
USA with lat=[20,50] and lon=[-130,-50]

Answers (1)

KSSV
KSSV on 9 Jan 2019
Read about interp2. This will help you to get the required data.
Let X,Y, A be your whole globe lat,lon and respective data.
lat=[20,50];
lon=[-130,-50] ;
M = 100 ; N = 100 ;
xi = linspace(lon(1),lon(2),M) ;
yi = linspace(lat(1),lat(2),N) ;
[Xi,Yi] = meshgrid(xi,yi) ;
Ai = interp2(X,Y,A,Xi,Yi) ;

Categories

Find more on Data Type Conversion 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!