fix problems with and duplication of MATLAB tools for EASEv2 latlon2i…

…nd (#640)

src/Applications/LDAS_App/util/inputs/mwRTM_params/get_mwRTM_vegcls_based.m

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

src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_ij_ind_from_latlon.m

@@ -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
@@ -21,4 +21,4 @@
 
 end
 
-%=========================================================================
+%=========================================================================

src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats.m

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

src/Applications/LDAS_App/util/postproc/write_smapL4SMqa.m

@@ -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:
@@ -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');
@@ -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 ) 

src/Applications/LDAS_App/util/shared/matlab/EASEv2_ind2latlon.m

@@ -198,3 +198,5 @@
     lat(idx) = NaN;
     lon(idx) = NaN;
 end
+
+% ========= EOF =========================================================

src/Applications/LDAS_App/util/shared/matlab/EASEv2_latlon2ind.m

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