How to extract data with latitude and longitude?

조회 수: 16 (최근 30일)
Trieu Hoa
Trieu Hoa 2019년 1월 9일
답변: KSSV 2019년 1월 9일
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]

답변 (1개)

KSSV
KSSV 2019년 1월 9일
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) ;

카테고리

Help CenterFile Exchange에서 Data Type Conversion에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by