/* * GVDR library for reading GVDR data files * Copyright (C) 1994 Michael J. Maurer * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Library General Public License for more details. * * You should have received a copy of the GNU Library General Public * License along with this library; if not, write to the Free * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. * * Michael Maurer * Durand Bldg - Room 232 * Stanford, CA 94305-4055 * (415) 723-1024 */ static char rcsid[]="$Id: callproj.c,v 1.2 1993/12/22 00:19:18 maurer Exp $"; /****************************************************************************** callproj.c Function: Direct interfact to USGS 'proj' cartographic library for converting between many cartesian map projections. Requires version 4 of 'proj'. For version 4.2.1 of PROJ, I changed line 15 of pj_fwd.c from if ((t = fabs(lp.phi)-HALFPI) > EPS || fabs(lp.lam) > 10.) { to if ((t = fabs(lp.phi)-HALFPI) > EPS || fabs(lp.lam) > 12.) { to give me more leeway in converting wraparound longitude. This file is part of the STARLab Magellan Altimeter Data Processing Software. Michael Maurer, November 1993. ******************************************************************************/ /* $Log: callproj.c,v $ * Revision 1.2 1993/12/22 00:19:18 maurer * Changed util.h to libmisc.h and added strdup() for mips. * * Revision 1.1 1993/11/19 01:32:18 maurer * Changed storage of ID number for PROJ interface. * * Revision 1.0 1993/11/03 21:04:20 maurer * Initial revision * */ #include #include #include #include "projects.h" /* proj library definitions */ #include "libmisc.h" #include "proj.h" #define CALLPROJ_C #include "callproj.p" /****************************************************************************** To localize knowledge of proj library, we keep all projections here in a local data structure. The outside world refers to them with an ID number indicating the position within the xform[] array. ******************************************************************************/ #define NXF 10 static int xfid = 0; static PJ *xform[NXF]; /****************************************************************************** callproj_init Initialize map projection transformation. The 'proj' library is initialized to convert from lon/lat pairs to x/y pairs (and vice versa), with details of the projection specified in P. Sets the field P->usgs_id to contain the transform ID number. Use callproj_query() to feed data to the process and received the transformed results. Returns nonzero on failure. ******************************************************************************/ int callproj_init(proj_t *P) { int argc,i; char *argv[20],buf[128]; P->usgs_id = -1; /* indicates failure */ if (xfid==0) for (i=0; i=NXF) return 1; /* build parameter list to initialize USGS proj */ argc=0; argv[argc++]="over"; /* allow "overflow" (no reduction to +-180) */ argv[argc++]="no_defs"; /* don't use definition files */ /* translate projection to name understood by USGS proj */ switch (P->proj) { case pjSinusoidal: argv[argc++]="proj=sinu"; break; case pjMercator: argv[argc++]="proj=merc"; break; case pjPolarStereo: case pjStereographic: argv[argc++]="proj=stere"; break; default: error(0,0,"[callproj_init] projection %d unsupported",P->proj); return 2; } /* only support square pixels at the moment */ if (P->par.ypixsiz!=P->par.xpixsiz) { error(0,0,"[callproj_init] ypixsiz must be same as xpixsiz"); return 3; } /* specify radius or major/minor axes of planet. scale these by the pixel size to duplicate "-m 1:xpixsiz" behavior. */ if (P->par.pradius < P->par.eradius) { sprintf(buf,"a=%f",P->par.eradius / P->par.xpixsiz); argv[argc++]=strdup(buf); sprintf(buf,"b=%f",P->par.pradius / P->par.xpixsiz); argv[argc++]=strdup(buf); } else if (P->par.pradius == P->par.eradius) { sprintf(buf,"R=%f",P->par.eradius / P->par.xpixsiz); argv[argc++]=strdup(buf); } else if (P->par.pradius > P->par.eradius) { error(0,0,"[callproj_init] can only handle oblate ellipsoids"); return 4; } /* add projection parameters */ sprintf(buf,"lon_0=%f",P->par.lon0); argv[argc++]=strdup(buf); sprintf(buf,"lat_0=%f",P->par.lat0); argv[argc++]=strdup(buf); sprintf(buf,"x_0=%f",P->par.x0); argv[argc++]=strdup(buf); sprintf(buf,"y_0=%f",P->par.y0); argv[argc++]=strdup(buf); argv[argc]=NULL; /* debuga(argv,%s,argc);*/ /* initialize projection */ if ((xform[xfid]=pj_init(argc,argv))==NULL) { error(0,0,"[callproj_init] pj_init failed \"%s\"",pj_strerrno(pj_errno)); return 5; } P->usgs_id=xfid; xfid++; return 0; } /****************************************************************************** callproj_query Feed data in xy[] to 'proj' library projection with ID number id. The data in xy[] is lon-lat or x-y in alternation, for a total of nxy pairs. On return, the data has been transformed in place. Units of latitude and longitude are radians; units of cartesian coordinates are pixel units. If inv!=0, the inverse transformation (that is, x-y to lon-lat) will be computed. Returns nonzero on failure. ******************************************************************************/ int callproj_query(int id, int inv, double *xy, int nxy) { UV data; int i; if (id<0 || id>=xfid || xform[id]==NULL) /* is this ID initialized? */ return 1; for (i=0; i=xfid || xform[id]==NULL) /* is this ID initialized? */ return 1; pj_free(xform[id]); xform[id]=NULL; return 0; } #ifdef mips char * strdup(s1) char *s1; { char *s2=NULL; if (s1) if (s2=malloc(strlen(s1)+1)) strcpy(s2,s1); return s2; } #endif