MtkLatLonToLSAry.c

Go to the documentation of this file.

/*===========================================================================
=                                                                           =
=                             MtkLatLonToLSAry                              =
=                                                                           =
=============================================================================

                         Jet Propulsion Laboratory
                                   MISR
                               MISR Toolkit

            Copyright 2005, California Institute of Technology.
                           ALL RIGHTS RESERVED.
                 U.S. Government Sponsorship acknowledged.

============================================================================*/

#include "MisrMapQuery.h"
#include "MisrCoordQuery.h"
#include "MisrUnitConv.h"
#include "MisrError.h"
#include "proj.h"
#include "gctp_prototypes.h"
#include "math.h"

MTKt_status MtkLatLonToLSAry(
  MTKt_MapInfo mapinfo, 
  int nelement,         
  const double *lat_dd, 
  const double *lon_dd, 
  float *line,          
  float *sample          )
{
  MTKt_status status_code;      /* Return status of this function */
  MTKt_status status;           /* Return status */
  int i;                        /* Loop index */
  MTKt_MisrProjParam pp;    /* Projection parameters */
  int iflg;         /* GCTP status flag */
  int (*for_trans[MAXPROJ+1])(double, double, double *, double *); /* Array of function ptrs (Not used) */
  double deg2rad = acos(-1) / 180.0;

  if (lat_dd == NULL || lon_dd == NULL || line == NULL || sample == NULL)
    MTK_ERR_CODE_JUMP(MTK_NULLPTR);

  if (nelement < 0)
    MTK_ERR_CODE_JUMP(MTK_BAD_ARGUMENT);

  status_code = MTK_SUCCESS;

  status = MtkPathToProjParam(mapinfo.path, 0, &pp);
  MTK_ERR_COND_JUMP(status);

  for_init(pp.projcode, pp.zonecode, pp.projparam, pp.spherecode,
       NULL, NULL, &iflg, for_trans);
  if (iflg)
    MTK_ERR_CODE_JUMP(MTK_GCTP_FORWARD_INIT_FAILED);

  for (i = 0; i < nelement; i++) {
    double som_x;                 /* SOM X */
    double som_y;                 /* SOM Y */
    float tline;
    float tsample; 

    iflg = somfor(lon_dd[i]*deg2rad, lat_dd[i]*deg2rad, &som_x, &som_y);
    if (iflg)
      MTK_ERR_CODE_JUMP(MTK_GCTP_FORWARD_PROJ_FAILED);

    tline = (float)((som_x - mapinfo.som.ulc.x) / mapinfo.resolution);
    tsample = (float)((som_y - mapinfo.som.ulc.y) / mapinfo.resolution);

    /* Check line/sample bounds */
    if (tline < -0.5 || tline > mapinfo.nline - 0.5) {
      tline = -1.0;
      tsample = -1.0;
      status_code = MTK_OUTBOUNDS;
    }

    if (tsample < -0.5 || tsample > mapinfo.nsample - 0.5) {
      tline = -1.0;
      tsample = -1.0;
      status_code = MTK_OUTBOUNDS;
    }

    line[i] = tline;
    sample[i] = tsample;
  }

  return status_code;

ERROR_HANDLE:
  return status_code;
}