/* * 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: proj.c,v 1.7 1994/04/29 02:02:00 maurer Exp $"; /****************************************************************************** proj.c Function: Provide basic map projection interface. This file is part of the STARLab Magellan Altimeter Data Processing Software. Michael Maurer, May 1993. ******************************************************************************/ /* $Log: proj.c,v $ * Revision 1.7 1994/04/29 02:02:00 maurer * Removed some unnecessary tests. * * Revision 1.6 1993/12/22 00:20:16 maurer * Split into proj.c and wproj.c. All write functions in wproj.c. * Changed util.h to libmisc.h for proj.c. * * Revision 1.5 1993/11/24 18:36:07 maurer * Moved gidx_clear() and gidx_free() to proj_clear() and proj_free(). * * Revision 1.4 1993/11/19 19:28:00 maurer * Bug fix: hangle possible null P->lin in proj_bounds(). * * Revision 1.3 1993/11/19 01:37:20 maurer * Now proj_pixlim() finds the maximum of the pixel extents for the * center, top and bottom of the line. Now any latitude or longitude * within the limits will pass the proj_bounds() test also. The * previous method could miss small wedges. * Now proj_alloc() can allocate memory for longitude/latitude values * of each pixel, and for additional user data. * Changed proj_getid() to store the USGS module's ID number inside the * proj_t structure, instead of having to search through all module * IDs and compare the projection parameters. * Added proj_lldist() and proj_xydist(). * * Revision 1.2 1993/11/03 21:06:02 maurer * Updated USGS interface to use PROJ version 4 (in module callproj.c), * instead of version 3 (in module pipeproj.c). * * Revision 1.1 1993/08/25 21:23:23 maurer * Added proj_pdsname(). * * Revision 1.0 1993/06/12 00:39:19 maurer * Initial revision * */ #include #include #include #include #include "libmisc.h" #include "proj.h" #define PROJ_C #include "proj.p" #include "callproj.p" /* Long and short names of the projections in proj_id enum */ static char *proj_names[pjNone+1][2] = { "Tiepoint", "tipt", "Gnomonic", "gnom", "Stereographic", "ster", "Orthographic", "orth", "Equidistant", "eqds", "Equivalent", "eqiv", "Lambert Conic1", "lcn1", "Lambert Conic2", "lcn2", "Mercator", "merc", "Polar Stereographic", "post", "Mollweide", "moll", "Albers1", "alb1", "Albers2", "alb2", "Lambert Cylindrical", "lcyl", "Lambert Polar", "lpol", "Bonne", "bonn", "Sinusoidal", "sinu", "Werner", "wern", "User-Defined", "user", "Aitoff", "aitf", "Hammer", "hamm", "Briesemeister", "brie", "Cylindrical", "cyln", "Unknown projection", "unkn" }; /* Names of the projections in hemi_id enum */ static char *hemi_names[] = { "Global", "Equatorial", "North", "South", "East", "West", "None" }; char * proj_name(proj_id proj) { if (proj>pjNone) proj=pjNone; return proj_names[proj][0]; } char * proj_sname(proj_id proj) { if (proj>pjNone) proj=pjNone; return proj_names[proj][1]; } char * hemi_name(hemi_id hemi) { if (hemi>hmNone) hemi=hmNone; return hemi_names[hemi]; } hemi_id hemi_lookup(char *name) { hemi_id hemi; char *strcasecmp(); if (name==NULL) return hmNone; for (hemi=0; hemiproj=pjNone; P->hemi=hmNone; P->linlo=0; P->linhi= -1; P->lin=NULL; P->usgs_id= -1; memset(&P->par,0,sizeof(P->par)); } /****************************************************************************** proj_free Frees the memory associated with map projection P. ******************************************************************************/ int proj_free(proj_t *P) { int i,j; if (P->lin) { for (i=P->linlo; i<=P->linhi; i++) { for (j=P->lin[i].pixlo; j<=P->lin[i].pixhi; j++) { if (P->lin[i].pix[j].Nal>0) free(P->lin[i].pix[j].f); P->lin[i].pix[j].Nf=0; P->lin[i].pix[j].Nal=0; P->lin[i].pix[j].f=NULL; } free(P->lin[i].pix + P->lin[i].pixlo); P->lin[i].pix=NULL; P->lin[i].line=0; P->lin[i].pixlo=0; P->lin[i].pixhi=0; } free(P->lin + P->linlo); } proj_clear(P); return 0; } /****************************************************************************** proj_bounds Tests if the point X/Y (in pixel units) lies within the bounds of the projection P. Returns 0 if the point is within bounds, 1 if the line Y is out of bounds, and 2 if the pixel X is out of bounds for the given line Y. ******************************************************************************/ int proj_bounds(proj_t *P, int X, int Y) { if (P->lin==NULL) return -1; if (Ylinlo || Y>P->linhi) return 1; if (Xlin[Y].pixlo || X>P->lin[Y].pixhi) return 2; return 0; } /****************************************************************************** proj_getid Looks up the desired cartographic projection in a static list. Returns the ID number of the projection assigned by callproj_init(). If an exact match is not found, callproj_init() is called to initialize the transformation software. ******************************************************************************/ static int proj_getid(proj_t *P) { int err; if (P->usgs_id<0) if (err=callproj_init(P)) error(-1,errno,"[proj_getid] callproj_init failed %d",err); return P->usgs_id; } /****************************************************************************** proj_fromlat Transforms a lon/lat pair (in degrees) into an X/Y pair (in pixel units) for the projection in P. Rather than perform computations here, the callproj.c module handles the interface to the USGS cartographic library 'proj'. The first time a particular projection is seen, the function proj_getid() sets up the communication channel and saves it in a static area. Subsequent calls communicate with the library and make transformation requests. ******************************************************************************/ int proj_fromlat(proj_t *P, double lon, double lat, double *X, double *Y) { int err,id; double xy[2]; id=proj_getid(P); #if 0 while (lon>360.0) /* allow both 0.0 and 360.0 */ lon-=360.0; while (lon<0.0) lon+=360.0; #endif xy[0]=lon*M_PI/180; xy[1]=lat*M_PI/180; if (err=callproj_query(id,0,xy,1)) { error(-1,errno,"[proj_fromlat] callproj_query failed %d",err); return 1; } *X=xy[0]; *Y=xy[1]; return 0; } /****************************************************************************** proj_tolat Transforms a X/Y pair (in pixel units) into an lon/lat pair (in degrees) for the projection in P. Rather than perform computations here, the callproj.c module handles the interface to the USGS cartographic library 'proj'. The first time a particular projection is seen, the function proj_getid() sets up the communication channel and saves it in a static area. Subsequent calls communicate with the library and make transformation requests. ******************************************************************************/ int proj_tolat(proj_t *P, double X, double Y, double *lon, double *lat) { int err,id,save=0; double uv[2]; /* shortcut if we've already computed and stored this */ if (X==(int)X && Y==(int)Y && proj_bounds(P,(int)X,(int)Y)==0 && P->lin[(int)Y].lon && P->lin[(int)Y].lat) { /* if we haven't computed it yet, save it after doing so */ save=1; if (P->lin[(int)Y].lon[(int)X] != HUGE_VAL && P->lin[(int)Y].lat[(int)X] != HUGE_VAL) { /* we have already computed it: just look it up */ *lon = P->lin[(int)Y].lon[(int)X]; *lat = P->lin[(int)Y].lat[(int)X]; return 0; } } id=proj_getid(P); uv[0]=X; uv[1]=Y; if (err=callproj_query(id,1,uv,1)) { error(-1,errno,"[proj_tolat] callproj_query failed %d",err); return 1; } *lon=uv[0]*180/M_PI; *lat=uv[1]*180/M_PI; /* if we are saving lon/lat, store the result now */ if (save) { P->lin[(int)Y].lon[(int)X] = *lon; P->lin[(int)Y].lat[(int)X] = *lat; } return 0; } /****************************************************************************** proj_phdr Print human readable dump of secondary projection header in P. ******************************************************************************/ void proj_phdr(FILE *fp, proj_t *P) { if (P->proj>pjNone) { fprintf(fp,"Unknown projection: %d",P->proj); return; } fprintf(fp," Projection: %s (%s)\n", proj_name(P->proj),hemi_name(P->hemi)); fprintf(fp,"X-Pixel Size: %8.1f X-offset: %+8.2f\n", P->par.xpixsiz,P->par.x0); fprintf(fp,"Y-Pixel Size: %8.1f Y-offset: %+8.2f\n", P->par.xpixsiz,P->par.y0); fprintf(fp," First Line: %8d E-Radius: %8.0f\n", P->linlo,P->par.eradius); fprintf(fp," Last Line: %8d P-Radius: %8.0f\n", P->linhi,P->par.pradius); fprintf(fp," Longitude_0: %+8.2f Latitude_0: %+8.2f\n", P->par.lon0,P->par.lat0); fprintf(fp," Longitude_1: %+8.2f Latitude_1: %+8.2f\n", P->par.lon1,P->par.lat1); fprintf(fp," Longitude_2: %+8.2f Latitude_2: %+8.2f\n", P->par.lon2,P->par.lat2); }