How to extract data with latitude and longitude?

Question

Trieu Hoa 2019-1-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/438864-how-to-extract-data-with-latitude-and-longitude

回答： 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.

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]

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

KSSV 2019-1-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/438864-how-to-extract-data-with-latitude-and-longitude#answer_355608

在 MATLAB Online 中打开

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) ;