Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fix problems and duplication with MATLAB tools for EASEv2 latlon2ind #640

Merged
merged 1 commit into from
Apr 12, 2023
Merged

fix problems and duplication with MATLAB tools for EASEv2 latlon2ind #640

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/Applications/LDAS_App/util/inputs/mwRTM_params/get_mwRTM_vegcls_based.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,11 +48,11 @@
%even if re-ordered output from an LDASsa-run is used as input, the
%resulting mwRTMparam.xxxx(:) is always monotonically increasingly sorted

[row_M09,col_M09] = smapeasev2_latlon2ind(tile_coord.com_lat,tile_coord.com_lon,'M09');
[row_M09,col_M09] = EASEv2_latlon2ind(tile_coord.com_lat,tile_coord.com_lon,'M09',1);
row_M09 = row_M09 +1;
col_M09 = col_M09 +1;

[row_M36,col_M36] = smapeasev2_latlon2ind(tile_coord.com_lat,tile_coord.com_lon,'M36');
[row_M36,col_M36] = EASEv2_latlon2ind(tile_coord.com_lat,tile_coord.com_lon,'M36',1);
row_M36 = row_M36 + 1;
col_M36 = col_M36 + 1;

Expand Down
6 changes: 3 additions & 3 deletions src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_ij_ind_from_latlon.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,11 @@
if (strcmp(tile_grid.gridtype,'EASEv2_M36'))
%row, col
[j_indg,i_indg] = ...
smapeasev2_latlon2ind(lat,lon,'M36');
EASEv2_latlon2ind(lat,lon,'M36',1);
elseif (strcmp(tile_grid.gridtype,'EASEv2_M09'))
%row, col
[j_indg,i_indg] = ...
smapeasev2_latlon2ind(lat,lon,'M09');
EASEv2_latlon2ind(lat,lon,'M09',1);
else
error('not ready for this grid');
end
Expand All @@ -21,4 +21,4 @@

end

%=========================================================================
%=========================================================================
5 changes: 3 additions & 2 deletions src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -231,9 +231,10 @@
%2) convert back to lat/lon at center of obs
if (~isempty(strfind(convert_grid, 'M36')) && ~isempty(strfind(convert_grid, 'EASEv2')))
gridid = 'M36';
[central_row,central_col] = smapeasev2_latlon2ind(central_lat,central_lon,gridid);
[central_lat,central_lon] = smapeasev2_ind2latlon(central_row,central_col,gridid);
[central_row,central_col] = EASEv2_latlon2ind(central_lat,central_lon,gridid,1);
[central_lat,central_lon] = EASEv2_ind2latlon(central_row,central_col,gridid);
elseif (~isempty(strfind(convert_grid, 'M36')) && ~isempty(strfind(convert_grid, 'EASEv1')))
error('Must provide smapeasev1_latlon2ind() and smapeasev1_ind2latlon()!')
gridid = 'M36';
[central_row,central_col] = smapeasev1_latlon2ind(central_lat,central_lon,gridid);
[central_lat,central_lon] = smapeasev1_ind2latlon(central_row,central_col,gridid);
Expand Down
223 changes: 5 additions & 218 deletions src/Applications/LDAS_App/util/postproc/write_smapL4SMqa.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,9 @@
function [] = write_smapL4SMqa( gph_aup_lmc_fnames, tilecoord_fname, tilegrids_fname )

% THE FOLLOWING PATH SHOULD BE ADDED IN MATLAB SCRIPT THAT CALLS THIS FUNCTION
% add path to matlab functions in src/Applications/LDAS_App/util/shared/matlab/
addpath('../shared/matlab/');

% Generate *.qa files from SMAP L4_SM "gph" or "aup" granules.
%
% Inputs:
Expand Down Expand Up @@ -452,7 +456,7 @@
% Map M09 tiles to M36 grid cells to verify the M36 obs

[M36_row,M36_col] = ...
smapeasev2_latlon2ind(tile_coord.com_lat,tile_coord.com_lon,'M36');
EASEv2_latlon2ind(tile_coord.com_lat,tile_coord.com_lon,'M36',1);

M36_row_col = [M36_row M36_col];
[unique_rc, ia, ic] = unique(M36_row_col,'rows');
Expand Down Expand Up @@ -1241,223 +1245,6 @@
array = array_out;
end

% *********************************************************************************************
% *********************************************************************************************
% *********************************************************************************************
%
% function [row,col] = smapeasev2_latlon2ind(lat,lon,gridid)
%
% SMAPEASE2FORWARD The principal function is to perform forward transformation
% from (lat,lon)'s to (row,col)'s for a set of nested EASE
% grids defined at 1, 3, 9, and 36km grid resolutions. These
% grids are all based on the EASE-Grid 2.0 specification (WGS84
% ellipsoid).
%
% SYNTAX [row,col] = smapease2forward(lat,lon,gridid)
%
% where gridid is a 3-character string enclosed in single
% quotes, in the form of {M|N|S}{01,03,09,36}. This subroutine
% accepts vector inputs and produce vector outputs.
%
% HISTORY This subroutine was adapted from the offical EASE-Grid-2.0
% conversion utilities (written in IDL) developed by the
% NSIDC.
%
% Note that in NSIDC's original implementation, (row,col) are
% zero-based. In other words, the first cell is (0,0) and the
% last cell is (N-1,M-1), where N and M are the row and column
% dimensions of the array. In this MATLAB implementation, the
% same convention is used. In other words, the end point of
% the first cell is located at (r,c) = (-0.5,-0.5) whereas the
% end point of the last cell is located at (r,c) = (14615.5,
% 34703.5). Thus,
%
% [lat,lon] = smapease2inverse(-0.5,-0.5,'M01') returns:
% lat = 85.044566407398861
% lon = 1.799999999999994e+02
%
% [lat,lon] = smapease2inverse(14615.5,34703.5,'M01') returns:
% lat = -85.044566407398861
% lon = -1.799999999999994e+02
%
% The polar grids, on the other hand, are more complete in
% terms of latitude coverage:
%
% [lat,lon] = smapease2inverse(8999,8999,'N01')
% lat = 89.993669248945238
% lon = -135
% [lat,lon] = smapease2inverse(9000,9000,'N01')
% lat = 89.993669248945238
% lon = 45
%
% [lat,lon] = smapease2inverse(8999,8999,'S01')
% lat = -89.993669248945238
% lon = -45
% [lat,lon] = smapease2inverse(9000,9000,'S01')
% lat = -89.993669248945238
% lon = 135
%
% UPDATE North/south polar projections were added. (03/2012)
%
% REFERENCE Brodzik, M. J., B. Billingsley, T. Haran, B. Raup, and M. H.
% Savoie (2012): EASE-Grid 2.0: Incremental but Significant
% Improvements for Earth-Gridded Data Sets. ISPRS International
% Journal of Geo-Information, vol. 1, no. 1, pp. 32-45,
% http://www.mdpi.com/2220-9964/1/1/32/
%
% Steven Chan, 11/2011
% Email: steven.k.chan@jpl.nasa.gov

function [row,col] = smapeasev2_latlon2ind(lat,lon,gridid)

% Constants returned by EASE2_GRID_INFO.PRO
projection = gridid(1);
switch gridid
case 'M36'
map_scale_m = 36032.220840584;
cols = 964;
rows = 406;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'M09'
map_scale_m = 9008.055210146;
cols = 3856;
rows = 1624;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'M03'
map_scale_m = 3002.6850700487;
cols = 11568;
rows = 4872;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'M01'
map_scale_m = 1000.89502334956;
cols = 34704;
rows = 14616;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'N36'
map_scale_m = 36000.0;
cols = 500;
rows = 500;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'N09'
map_scale_m = 9000.0;
cols = 2000;
rows = 2000;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'N03'
map_scale_m = 3000.0;
cols = 6000;
rows = 6000;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'N01'
map_scale_m = 1000.0;
cols = 18000;
rows = 18000;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'S36'
map_scale_m = 36000.0;
cols = 500;
rows = 500;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'S09'
map_scale_m = 9000.0;
cols = 2000;
rows = 2000;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'S03'
map_scale_m = 3000.0;
cols = 6000;
rows = 6000;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
case 'S01'
map_scale_m = 1000.0;
cols = 18000;
rows = 18000;
r0 = (cols-1)/2;
s0 = (rows-1)/2;
otherwise
disp(['ERROR: Incompatible grid specification.']);
end

% Constants returned by EASE2_MAP_INFO.PRO
epsilon = 1.0e-6;
map_equatorial_radius_m = 6378137.0;
map_eccentricity = 0.081819190843;
e2 = map_eccentricity^2;
switch projection
case 'M'
map_reference_latitude = 0.0;
map_reference_longitude = 0.0;
map_second_reference_latitude = 30.0;
sin_phi1 = sin(map_second_reference_latitude*pi/180);
cos_phi1 = cos(map_second_reference_latitude*pi/180);
kz = cos_phi1/sqrt(1.0-e2*sin_phi1*sin_phi1);
case 'N'
map_reference_latitude = 90.0;
map_reference_longitude = 0.0;
case 'S'
map_reference_latitude = -90.0;
map_reference_longitude = 0.0;
end

% Selected calculations inside WGS84_CONVERT.PRO and WGS84_CONVERT_XY.PRO
dlon = lon - map_reference_longitude;
msk1 = dlon < -180.0; dlon(msk1) = dlon(msk1) + 360.0;
msk2 = dlon > +180.0; dlon(msk2) = dlon(msk2) - 360.0;
phi = lat*pi/180.0;
lam = dlon*pi/180.0;
sin_phi = sin(phi);
q = (1.0-e2)*((sin_phi./(1.0-e2*sin_phi.*sin_phi))-(1.0/(2.0*map_eccentricity))*log((1.0-map_eccentricity*sin_phi)./(1.0+map_eccentricity*sin_phi)));
qp = 1.0-((1.0-e2)/(2.0*map_eccentricity)*log((1.0-map_eccentricity)/(1.0+map_eccentricity)));
switch projection
case 'M'
x = map_equatorial_radius_m*kz*lam;
y = (map_equatorial_radius_m*q)/(2.0*kz);
case 'N'
tmp = qp - q;
tmp(abs(tmp) < epsilon) = 0.0;
rho = map_equatorial_radius_m*sqrt(tmp);
x = rho.*sin(lam);
y = -rho.*cos(lam);
case 'S'
tmp = qp + q;
tmp(abs(tmp) < epsilon) = 0.0;
rho = map_equatorial_radius_m*sqrt(tmp);
x = rho.*sin(lam);
y = rho.*cos(lam);
end
row = s0-(y/map_scale_m);
col = r0+(x/map_scale_m);
switch projection
case 'N'
idx = lat < 0.0;
row(idx) = NaN;
col(idx) = NaN;
case 'S'
idx = lat > 0.0;
row(idx) = NaN;
col(idx) = NaN;
end

% restrict indices to (zero-based) min and max indices
% (needed because of round-off error - GDL, 07 Jun 2013)

row = max( 0, min( rows-1, round(row)));
col = max( 0, min( cols-1, round(col)));


% *********************************************************************************************
% *********************************************************************************************
% *********************************************************************************************

function [ stats_out ] = stats_array( array_in, weights, Nstat )
Expand Down
2 changes: 2 additions & 0 deletions src/Applications/LDAS_App/util/shared/matlab/EASEv2_ind2latlon.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,3 +198,5 @@
lat(idx) = NaN;
lon(idx) = NaN;
end

% ========= EOF =========================================================
18 changes: 17 additions & 1 deletion src/Applications/LDAS_App/util/shared/matlab/EASEv2_latlon2ind.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,14 @@
% Steven Chan, 11/2011
% Email: steven.k.chan@jpl.nasa.gov

function [row,col] = EASEv2_latlon2ind(lat,lon,gridid)
function [row,col] = EASEv2_latlon2ind(lat,lon,gridid,return_rounded)

% By design, [row, col] are real numbers, with the fractional portion indicating
% the position of the specified [lat, lon] coordinates between adjacent grid
% cell centers. E.g., col=0.5 indicates that the input longitude is on the
% boundary between grid cells associated with col=0 and col=1.
% If the [optional] input argument 'return_rounded' is present and ~=0, then
% [row, col] are rounded to the nearest integer.

% Constants returned by EASE2_GRID_INFO.PRO
projection = gridid(1);
Expand Down Expand Up @@ -200,3 +207,12 @@
row(idx) = NaN;
col(idx) = NaN;
end

if exist('return_rounded','var')
if return_rounded
col=round(col);
row=round(row);
end
end

% ========= EOF =========================================================