PLplot  5.10.0
plmap.c
Go to the documentation of this file.
00001 //      Continental Outline and Political Boundary Backgrounds
00002 //
00003 //      Some plots need a geographical background such as the global
00004 //      surface temperatures or the population density.  The routine
00005 //      plmap() will draw one of the following backgrounds: continental
00006 //      outlines, political boundaries, the United States, and the United
00007 //      States with the continental outlines.  The routine plmeridians()
00008 //      will add the latitudes and longitudes to the background.  After
00009 //      the background has been drawn, one can use a contour routine or a
00010 //      symbol plotter to finish off the plot.
00011 //
00012 //      Copyright (C) 1991, 1993, 1994  Wesley Ebisuzaki
00013 //      Copyright (C) 1994, 2000, 2001  Maurice LeBrun
00014 //      Copyright (C) 1999  Geoffrey Furnish
00015 //      Copyright (C) 2000-2014 Alan W. Irwin
00016 //      Copyright (C) 2001  Andrew Roach
00017 //      Copyright (C) 2001, 2004  Rafael Laboissiere
00018 //      Copyright (C) 2002  Vincent Darley
00019 //      Copyright (C) 2003  Joao Cardoso
00020 //
00021 //      This file is part of PLplot.
00022 //
00023 //      PLplot is free software; you can redistribute it and/or modify
00024 //      it under the terms of the GNU Library General Public License
00025 //      as published by the Free Software Foundation; version 2 of the
00026 //      License.
00027 //
00028 //      PLplot is distributed in the hope that it will be useful, but
00029 //      WITHOUT ANY WARRANTY; without even the implied warranty of
00030 //      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00031 //      GNU Library General Public License for more details.
00032 //
00033 //      You should have received a copy of the GNU Library General Public
00034 //      License along with this library; if not, write to the Free Software
00035 //      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301
00036 //      USA
00037 //
00038 //
00039 
00040 #define DEBUG
00041 #define NEED_PLDEBUG
00042 
00043 #include "plplotP.h"
00044 
00045 #ifdef HAVE_SHAPELIB
00046 #include <shapefil.h>
00047 
00048 SHPHandle
00049 OpenShapeFile( const char *fn );
00050 
00051 #ifdef HAVE_SAHOOKS
00052 static void
00053 CustomErrors( const char *message );
00054 #endif //HAVE_SAHOOKS
00055 
00056 #define MAP_FILE    ""
00057 #define OpenMap     OpenShapeFile
00058 #define CloseMap    SHPClose
00059 
00060 #else //HAVE_SHAPELIB
00061 
00062 #define MAP_FILE    ".map"
00063 #define OpenMap     plLibOpenPdfstrm
00064 #define CloseMap    pdf_close
00065 #define OFFSET      ( 180 * 100 )
00066 #define SCALE       100.0
00067 #define W_BUFSIZ    ( 32 * 1024 )
00068 #define SHPT_ARC     1
00069 #define SHPT_POINT   2
00070 #define SHPT_POLYGON 3
00071 
00072 #endif //HAVE_SHAPELIB
00073 
00074 //redistributes the lon value onto either 0-360 or -180-180 for wrapping
00075 //purposes.
00076 void
00077 rebaselon( PLFLT *lon, PLFLT midlon )
00078 {
00079         if( *lon > midlon + 180.0 )
00080                 *lon -= floor( ( *lon - midlon - 180.0 ) / 360.0 + 1.0 ) * 360.0;
00081         else if( *lon < midlon - 180.0 )
00082                 *lon += floor( ( midlon - 180.0 - *lon ) / 360.0 + 1.0 ) * 360.0;
00083         
00084 }
00085 
00086 //append a PLFLT to an array of PLFLTs. array is a pointer to the array,
00087 //n is the current size of the array, val is the value to append.
00088 //returns 0 for success, 1 for could not allocate new memory. If memory
00089 //could not be allocated then the previous array remains intact
00090 int appendflt( PLFLT **array, size_t n, PLFLT val )
00091 {
00092         size_t i;
00093         PLFLT *temp = ( PLFLT* ) malloc( ( n + 1 ) * sizeof( PLFLT ) );
00094         if( !temp)
00095                 return 1;
00096         for( i = 0; i < n; ++i )
00097                 temp[i]=(*array)[i];
00098         temp[n] = val;
00099         free ( *array );
00100         *array = temp;
00101         return 0;
00102 }
00103 
00104 //As for appendflt, but for an array of ints
00105 int appendint( int **array, size_t n, int val )
00106 {
00107         size_t i;
00108         int *temp = ( int* ) malloc( ( n + 1 ) * sizeof( int ) );
00109         if( !temp)
00110                 return 1;
00111         for( i = 0; i < n; ++i )
00112                 temp[i]=(*array)[i];
00113         temp[n] = val;
00114         free ( *array );
00115         *array = temp;
00116         return 0;
00117 }
00118 
00119 //As for appendflt, but for an array of pointers to PLFLTs
00120 int appendfltptr( PLFLT ***array, size_t n, PLFLT *val )
00121 {
00122         size_t i;
00123         PLFLT **temp = ( PLFLT** ) malloc( ( n + 1 ) * sizeof( PLFLT *) );
00124         if( !temp)
00125                 return 1;
00126         for(i = 0; i < n; ++i)
00127                 temp[i]=(*array)[i];
00128         temp[n] = val;
00129         free ( *array );
00130         *array = temp;
00131         return 0;
00132 }
00133 
00134 //Check to see if the mapform wraps round longitudes. For example, a polar
00135 // projection wraps round longitudes, but a cylindrical projection does not.
00136 //Returns 1 if the mapform wraps or 0 if not.
00137 char
00138 checkwrap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), PLFLT lon, PLFLT lat )
00139 {
00140         PLFLT x[]={lon};
00141         PLFLT y[]={lat};
00142         PLFLT resultx;
00143         PLFLT resulty;
00144 
00145         if( !mapform )
00146                 return 0;
00147 
00148         mapform( 1, x, y );
00149         resultx = x[0];
00150         resulty = y[0];
00151         x[0] = lon + 360;
00152         y[0] = lat;
00153         return ( ( ABS( x[0] - resultx ) < 1.0e-12 ) && ( ABS( y[0] - resulty ) < 1.0e-12 ) );
00154 }
00155 
00156 //-------------------------------------------------------------------------------------
00157 //Actually draw the map lines points and text.
00158 //-------------------------------------------------------------------------------------
00159 void
00160 drawmapdata( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), int shapetype, PLINT n, PLFLT *x, PLFLT *y, PLFLT dx, PLFLT dy, PLFLT just, const char *text )
00161 {
00162         PLINT i;
00163 
00164         //do the transform if needed
00165         if ( mapform != NULL )
00166                 ( *mapform )( n, x, y );
00167         
00168         if ( shapetype == SHPT_ARC )
00169                 plline( n, x, y );
00170         else if ( shapetype == SHPT_POINT )
00171                 for ( i = 0; i < n; ++i )
00172                         plptex( x[i], y[i], dx, dy, just, text );
00173         else if ( shapetype == SHPT_POLYGON )
00174                 plfill( n, x, y );
00175 #ifdef HAVE_SHAPELIB
00176         else if( shapetype == SHPT_ARCZ || shapetype == SHPT_ARCM )
00177                 plline( n, x, y );
00178         else if( shapetype == SHPT_POLYGON || shapetype == SHPT_POLYGONZ || shapetype == SHPT_POLYGONM )
00179                 plfill( n, x, y );
00180         else if( shapetype == SHPT_POINT || shapetype ==SHPT_POINTM || shapetype == SHPT_POINTZ )
00181                 for ( i = 0; i < n; ++i )
00182                         plptex( x[i], y[i], dx, dy, just, text );
00183 #endif //HAVE_SHAPELIB
00184 }
00185 
00186 
00187 //-------------------------------------------------------------------------------------
00188 //This is a function called by the front end map functions to do the map drawing. Its
00189 //parameters are:
00190 //mapform: The transform used to convert the data in raw coordinates to x, y positions
00191 //on the plot
00192 //type: either one of the plplot provided lat/lon maps or the path/file name of a 
00193 //shapefile
00194 //dx/dy: the gradient of text/symbols drawn if text is non-null
00195 //shapetype: one of ARC, SHPT_ARCZ, SHPT_ARCM, SHPT_POLYGON, SHPT_POLYGONZ, 
00196 //SHPT_POLYGONM, SHPT_POINT, SHPT_POINTM, SHPT_POINTZ. See drawmapdata() for the
00197 //how each type is rendered. But Basically the ARC options are lines, the POLYGON
00198 //options are filled polygons, the POINT options are points/text. Options beginning
00199 //SHPT will only be defined if HAVE_SHAPELIB is true
00200 //text: The text (which can be actual text or a unicode symbol) to be drawn at
00201 //each point
00202 //minlong/maxlong: The min/max longitude when using a plplot provided map or x value if
00203 //using a shapefile
00204 //minlat/maxlat: The min/max latitude when using a plplot provided map or y value if
00205 //using a shapefile
00206 //plotentries: used only for shapefiles, as one shapefile contains multiple vectors
00207 //each representing a different item (e.g. multiple boundaries, multiple height
00208 //contours etc. plotentries is an array containing the indices of the
00209 //entries within the shapefile that you wish to plot. if plotentries is null all 
00210 //entries are plotted
00211 //nplotentries: the number of elements in plotentries. Ignored if plplot was not built
00212 //with shapefile support or if plotentries is null
00213 //-------------------------------------------------------------------------------------
00214 void
00215 drawmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type,
00216        PLFLT dx, PLFLT dy, int shapetype, PLFLT just, const char *text,
00217        PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, int* plotentries, int nplotentries )
00218 {
00219 #if defined ( HAVE_SHAPELIB ) || defined ( PL_DEPRECATED )
00220     int       i, j;
00221     char      *filename = NULL;
00222     char      truncatedfilename[900];
00223     char      warning[1024];
00224     int       nVertices = 200;
00225     PLFLT     minsectlon, maxsectlon, minsectlat, maxsectlat;
00226      PLFLT     *bufx  = NULL, *bufy = NULL;
00227     int       bufsize = 0;
00228         int       filenamelen;
00229         PLFLT     **splitx = NULL;
00230         PLFLT     **splity = NULL;
00231         int       *splitsectionlengths = NULL;
00232         int       nsplitsections;
00233         PLFLT     lastsplitpointx;
00234         PLFLT     lastsplitpointy;
00235         PLFLT     penultimatesplitpointx;
00236         PLFLT     penultimatesplitpointy;
00237         char      islatlon = 1;
00238 
00239 
00240 #ifdef HAVE_SHAPELIB
00241     SHPHandle        in;
00242     int              nentries;
00243         int              entryindex=0;
00244     // Unnecessarily set nparts to quiet -O3 -Wuninitialized warnings.
00245     //int              nparts      = 0;
00246     int              entrynumber = 0;
00247     int              partnumber  = 0;
00248     double           mins[4];
00249     double           maxs[4];
00250     SHPObject        *object = NULL;
00251     double           *bufxraw;
00252     double           *bufyraw;
00253         char             *prjfilename = NULL;
00254         PDFstrm          *prjfile;
00255         char             prjtype[]={0, 0, 0, 0, 0, 0, 0};
00256         int              appendresult=0;
00257 #else
00258     PDFstrm          *in;
00259     //PLFLT            bufx[ncopies][200], bufy[ncopies][200];
00260     unsigned char    n_buff[2], buff[800];
00261     long int         t;
00262 #endif
00263 
00264     //
00265     // read map outline
00266     //
00267 
00268         //strip the .shp extension if a shapefile has been provided and add
00269         //the needed map file extension if we are not using shapefile
00270         if( strstr( type, ".shp" ) )
00271                 filenamelen=( int )( type-strstr( type, ".shp" ) );
00272         else
00273                 filenamelen = ( int )strlen( type );
00274     filename = ( char* )malloc( filenamelen + strlen( MAP_FILE ) + 1 );
00275         if(!filename)
00276         {
00277                 plabort("Could not allocate memory for concatenating map filename");
00278                 return;
00279         }
00280     strncpy( filename, type, filenamelen );
00281         filename[ filenamelen ] = '\0';
00282     strcat( filename, MAP_FILE );
00283 
00284         //copy the filename to a fixed length array in case it is needed for warning messages
00285         if( strlen( filename ) < 899 )
00286                 strcpy( truncatedfilename, filename );
00287         else
00288         {
00289                 memcpy( truncatedfilename, filename, 896 );
00290                 truncatedfilename[896] = '.';
00291                 truncatedfilename[897] = '.';
00292                 truncatedfilename[898] = '.';
00293                 truncatedfilename[899] = '\0';
00294         }
00295 
00296     strcpy( warning, "Could not find " );
00297     strcat( warning, filename );
00298     strcat( warning, " file." );
00299 #ifdef HAVE_SHAPELIB
00300         //Open the shp and shx file using shapelib
00301     if ( ( in = OpenShapeFile( filename ) ) == NULL )
00302     {
00303         plabort( warning );
00304                 free(filename);
00305         return;
00306     }
00307     SHPGetInfo( in, &nentries, &shapetype, mins, maxs );
00308         //also check for a prj file which will tell us if the data is lat/lon or projected
00309         //if it is projected then set ncopies to 1 - i.e. don't wrap round longitudes
00310     prjfilename = ( char* )malloc( filenamelen + 5 );
00311         if(!prjfilename)
00312         {
00313                 plabort("Could not allocate memory for generating map projection filename");
00314                 free(filename);
00315                 return;
00316         }
00317     strncpy( prjfilename, type, filenamelen );
00318         prjfilename[ filenamelen ] = '\0';
00319     strcat( prjfilename, ".prj" );
00320         prjfile = plLibOpenPdfstrm( prjfilename );
00321         if ( prjfile && prjfile->file )
00322         {
00323                 fread( prjtype, 1, 6, prjfile->file );
00324                 if( strcmp ( prjtype, "PROJCS" ) == 0 )
00325                         islatlon = 0;
00326                 pdf_close( prjfile );
00327         }
00328         free ( prjfilename );
00329         prjfilename = NULL;
00330 #else
00331     if ( ( in = plLibOpenPdfstrm( filename ) ) == NULL )
00332     {
00333         plwarn( warning );
00334         return;
00335     }
00336 #endif
00337 
00338     bufx = NULL;
00339     bufy = NULL;
00340 
00341     for (;; )
00342     {
00343 #ifdef HAVE_SHAPELIB
00344                 //each object in the shapefile is split into parts.
00345                 //If we are need to plot the first part of an object then read in a new object
00346                 //and check how many parts it has. Otherwise use the object->panPartStart vector
00347                 //to check the offset of this part and the next part and allocate memory. Copy
00348                 //the data to this memory converting it to PLFLT and draw it.
00349                 //finally increment the part number or if we have finished with the object reset the
00350                 //part numberand increment the object.
00351 
00352         //break condition if we've reached the end of the file
00353                 if ( ( !plotentries && ( entrynumber == nentries ) ) || ( plotentries && ( entryindex == nplotentries ) ) )
00354             break;
00355 
00356         //if partnumber == 0 then we need to load the next object
00357         if ( partnumber == 0 )
00358         {
00359                         if( plotentries )
00360                                 object = SHPReadObject( in, plotentries[entryindex] );
00361                         else
00362                                 object = SHPReadObject( in, entrynumber );
00363         }
00364                 //if the object could not be read, increment the object index to read and
00365                 //return to the top of the loop to try the next object.
00366                 if ( object == NULL )
00367         {
00368             entrynumber++;
00369                         entryindex++;
00370             partnumber = 0;
00371                         continue;
00372         }
00373 
00374         //work out how many points are in the current part
00375                 if( object->nParts == 0 )
00376                         nVertices = object->nVertices; //if object->nParts==0, we can still have 1 vertex. A bit odd but it's the way it goes
00377         else if ( partnumber == ( object->nParts - 1 ) )
00378             nVertices = object->nVertices - object->panPartStart[partnumber];//panPartStart holds the offset for each part
00379         else
00380             nVertices = object->panPartStart[partnumber + 1] - object->panPartStart[partnumber];//panPartStart holds the offset for each part
00381 #endif
00382         //allocate memory for the data
00383         if ( nVertices > bufsize )
00384         {
00385             bufsize = nVertices;
00386             free( bufx );
00387             free( bufy );
00388                         bufx = (PLFLT*)malloc( (size_t) bufsize * sizeof ( PLFLT ) );
00389                         bufy = (PLFLT*)malloc( (size_t) bufsize * sizeof ( PLFLT ) );
00390                         if(!bufx || ! bufy)
00391                         {
00392                                 plabort("Could not allocate memory for map data");
00393                                 free(filename);
00394                                 free (bufx);
00395                                 free (bufy);
00396                                 return;
00397                         }
00398         }
00399 
00400 #ifdef HAVE_SHAPELIB
00401         //point the plot buffer to the correct starting vertex
00402         //and copy it to the PLFLT arrays. If we had object->nParts == 0
00403                 //then panPartStart will be NULL
00404                 if(object->nParts>0)
00405                 {
00406                         bufxraw = object->padfX + object->panPartStart[partnumber];
00407                         bufyraw = object->padfY + object->panPartStart[partnumber];
00408                 }
00409                 else
00410                 {
00411                         bufxraw = object->padfX;
00412                         bufyraw = object->padfY;
00413                 }
00414 
00415         for ( i = 0; i < nVertices; i++ )
00416         {
00417             bufx[i] = (PLFLT) bufxraw[i];
00418             bufy[i] = (PLFLT) bufyraw[i];
00419         }
00420 
00421         //set the min x/y of the object
00422         minsectlon = object->dfXMin;
00423         maxsectlon = object->dfXMax;
00424         minsectlat = object->dfYMin;
00425         maxsectlat = object->dfYMax;
00426 
00427         //increment the partnumber or if we've reached the end of
00428         //an entry increment the entrynumber and set partnumber to 0
00429         if ( partnumber == object->nParts - 1 || object->nParts == 0 )
00430         {
00431             entrynumber++;
00432                         entryindex++;
00433             partnumber = 0;
00434                         SHPDestroyObject( object );
00435                         object=NULL;
00436         }
00437         else
00438             partnumber++;
00439 
00440                 if (nVertices==0)
00441                         continue;
00442 #else
00443         // read in # points in segment
00444         if ( pdf_rdx( n_buff, (long) sizeof ( unsigned char ) * 2, in ) == 0 )
00445             break;
00446         nVertices = ( n_buff[0] << 8 ) + n_buff[1];
00447         if ( nVertices == 0 )
00448             break;
00449 
00450         pdf_rdx( buff, (long) sizeof ( unsigned char ) * 4 * nVertices, in );
00451         if ( nVertices == 1 )
00452             continue;
00453 
00454         for ( j = i = 0; i < nVertices; i++, j += 2 )
00455         {
00456             t            = ( buff[j] << 8 ) + buff[j + 1];
00457             bufx[i] = ( (PLFLT) t - OFFSET ) / SCALE;
00458         }
00459         for ( i = 0; i < nVertices; i++, j += 2 )
00460         {
00461             t          = ( buff[j] << 8 ) + buff[j + 1];
00462             bufy[i] = ( (PLFLT) t - OFFSET ) / SCALE;
00463         }
00464         //set the min/max section lat/lon with extreme values
00465         //to be overwritten later
00466         minsectlon = 1000.;
00467         maxsectlon = -1000.;
00468         minsectlat = 1000.;
00469         maxsectlat = -1000.;
00470 
00471 #endif
00472 
00473                 if( islatlon )
00474                 {
00475                         //two obvious issues exist here with plotting longitudes:
00476                         //
00477                         //1) wraparound causing lines which go the wrong way round
00478                         //   the globe
00479                         //2) some people plot lon from 0-360 deg, others from -180 - +180
00480                         //
00481                         //we can cure these problems by conditionally adding/subtracting
00482                         //360 degrees to each data point in order to ensure that the
00483                         //distance between adgacent points is always less than 180
00484                         //degrees, then plotting up to 2 out of 5 copies of the data
00485                         //each separated by 360 degrees.
00486                 
00487                         //arrays of pointers to the starts of each section of data that
00488                         //has been split due to longitude wrapping, and an array of ints 
00489                         //to hold their lengths. Start with splitx and splity having one
00490                         //element pointing to the beginning of bufx and bufy
00491                         splitx = ( PLFLT** ) malloc ( sizeof ( PLFLT* ) );
00492                         splity = ( PLFLT** ) malloc ( sizeof ( PLFLT* ) );
00493                         //lengths of the split sections
00494                         splitsectionlengths = ( int* ) malloc ( sizeof ( size_t ) );
00495                         if(!splitx || ! splity || ! splitsectionlengths)
00496                         {
00497                                 plabort("Could not allocate memory for longitudinally split map data");
00498                                 free( filename );
00499                                 free ( bufx );
00500                                 free ( bufy );
00501                                 free ( splitx );
00502                                 free ( splity );
00503                                 free ( splitsectionlengths );
00504                                 return;
00505                         }
00506                         splitsectionlengths[0]=nVertices;
00507                         nsplitsections=1;
00508                         splitx[0]=bufx;
00509                         splity[0]=bufy;
00510 
00511                         //set the min/max lats/lons
00512                         minsectlon = MIN( minsectlon, bufx[0] );
00513                         maxsectlon = MAX( minsectlon, bufx[0] );
00514                         minsectlat = MIN( minsectlat, bufy[0] );
00515                         maxsectlat = MAX( maxsectlat, bufy[0] );
00516                         //ensure our lat and lon are on 0-360 grid and split the
00517                         //data where it wraps.
00518                         rebaselon( &bufx[0], ( minlong + maxlong ) / 2.0 );
00519                         for ( i = 1; i < nVertices; i++ )
00520                         {
00521                                 //put lon into 0-360 degree range
00522                                 rebaselon( &bufx[i], ( minlong + maxlong ) / 2.0 );
00523 
00524                                 //check if the previous point is more than 180 degrees away
00525                                 if ( bufx[i - 1] - bufx[i] > 180. || bufx[i - 1] - bufx[i] < -180.)
00526                                 {
00527                                         //check if the map transform deals with wrapping itself, e.g. in a polar projection
00528                                         //in this case give one point overlap to the sections so that lines are contiguous
00529                                         if( checkwrap( mapform, bufx[i], bufy[i] ) )
00530                                         {
00531                                                 appendresult += appendfltptr( &splitx, nsplitsections, bufx+i );
00532                                                 appendresult += appendfltptr( &splity, nsplitsections, bufy+i );
00533                                                 appendresult += appendint( &splitsectionlengths, nsplitsections, nVertices - i );
00534                                                 splitsectionlengths[nsplitsections-1] -= splitsectionlengths[nsplitsections] - 1;
00535                                                 nsplitsections++;
00536                                         }
00537                                         //if the transform doesn't deal with wrapping then allow 2 points overlap to fill in the 
00538                                         //edges
00539                                         else
00540                                         {
00541                                                 appendresult += appendfltptr( &splitx, nsplitsections, bufx+i-1 );
00542                                                 appendresult += appendfltptr( &splity, nsplitsections, bufy+i-1 );
00543                                                 appendresult += appendint( &splitsectionlengths, nsplitsections, nVertices - i + 1 );
00544                                                 splitsectionlengths[nsplitsections-1] -= splitsectionlengths[nsplitsections] - 2;
00545                                                 nsplitsections++;
00546                                         }
00547                                         if( appendresult > 0 )
00548                                         {
00549                                                 plabort("Could not allocate memory for appending to longitudinally split map data");
00550                                                 free( filename );
00551                                                 free ( bufx );
00552                                                 free ( bufy );
00553                                                 free ( splitx );
00554                                                 free ( splity );
00555                                                 free ( splitsectionlengths );
00556                                                 return;
00557                                         }
00558 
00559                                 }
00560 
00561                                 //update the mins and maxs
00562                                 minsectlon = MIN( minsectlon, bufx[i] );
00563                                 maxsectlon = MAX( minsectlon, bufx[i] );
00564                                 minsectlat = MIN( minsectlat, bufy[i] );
00565                                 maxsectlat = MAX( maxsectlat, bufy[i] );
00566                         }
00567 
00568                         //check if the latitude and longitude range means we need to plot this section
00569                         if ( ( maxsectlat > minlat ) && ( minsectlat < maxlat ) 
00570                                 && ( maxsectlon > minlong ) && ( minsectlon < maxlong ))
00571                         {
00572                                 //plot each split in turn, now is where we deal with the end points to
00573                                 //ensure we draw to the edge of the map
00574                                 for( i = 0; i < nsplitsections; ++i )
00575                                 {
00576                                         //check if the first 2 or last 1 points of the split section need 
00577                                         //wrapping and add or subtract 360 from them. Note that when the next
00578                                         //section is drawn the code below will undo this if needed
00579                                         if( splitsectionlengths[i] > 2 )
00580                                         {
00581                                                 if ( splitx[i][1] - splitx[i][2] > 180. )
00582                                                         splitx[i][1] -= 360.0;
00583                                                 else if ( splitx[i][1] - splitx[i][2] < -180. )
00584                                                         splitx[i][1] += 360.0;
00585                                         }
00586 
00587                                         if ( splitx[i][0] - splitx[i][1] > 180. )
00588                                                 splitx[i][0] -= 360.0;
00589                                         else if ( splitx[i][0] - splitx[i][1] < -180. )
00590                                                 splitx[i][0] += 360.0;
00591                                 
00592                                         if ( splitx[i][splitsectionlengths[i] - 2] - splitx[i][splitsectionlengths[i] - 1] > 180. )
00593                                                 splitx[i][splitsectionlengths[i] - 1] += 360.0;
00594                                         else if ( splitx[i][splitsectionlengths[i] - 2] - splitx[i][splitsectionlengths[i] - 1] < -180. )
00595                                                 splitx[i][splitsectionlengths[i] - 1] -= 360.0;
00596 
00597                                         //save the last 2 points - they will be needed by the next
00598                                         //split section and will be overwritten by the mapform
00599                                         lastsplitpointx = splitx[i][splitsectionlengths[i]-1];
00600                                         lastsplitpointy = splity[i][splitsectionlengths[i]-1];
00601                                         penultimatesplitpointx = splitx[i][splitsectionlengths[i]-2];
00602                                         penultimatesplitpointy = splity[i][splitsectionlengths[i]-2];
00603 
00604                                         //draw the split section
00605                                         drawmapdata(mapform, shapetype, splitsectionlengths[i], splitx[i], splity[i], dx, dy, just, text );
00606 
00607                                         for(j=1; j<splitsectionlengths[i]; ++j)
00608                                         {
00609                                                 if((splitx[i][j]<200.0 && splitx[i][j-1] > 260.0) || ( splitx[i][j-1]<200.0 && splitx[i][j] > 260.0 ))
00610                                                         plwarn("wrapping error");
00611                                         }
00612 
00613                                         //restore the last 2 points
00614                                         splitx[i][splitsectionlengths[i]-1] = lastsplitpointx;
00615                                         splity[i][splitsectionlengths[i]-1] = lastsplitpointy;
00616                                         splitx[i][splitsectionlengths[i]-2] = penultimatesplitpointx;
00617                                         splity[i][splitsectionlengths[i]-2] = penultimatesplitpointy;
00618                                 }
00619                         }
00620                 }
00621                 else
00622                 {
00623                         drawmapdata(mapform, shapetype, nVertices, bufx, bufy, dx, dy, just, text );
00624                 }
00625                 
00626                 free( splitx );
00627                 free( splity );
00628                 free( splitsectionlengths );
00629 
00630 
00631     }
00632     // Close map file
00633 #ifdef HAVE_SHAPELIB
00634     SHPClose( in );
00635 #else
00636     pdf_close( in );
00637 #endif
00638 
00639     //free memory
00640     free( bufx );
00641     free( bufy );
00642     free( filename );
00643 #else   // defined (HAVE_SHAPELIB) || defined (PL_DEPRECATED)
00644     plwarn( "Use of the old plplot map file format is deprecated.\nIt is recommended that the shapelib library be used to provide map support.\n" );
00645 #endif  // defined (HAVE_SHAPELIB) || defined (PL_DEPRECATED)
00646 }
00647 
00648 
00649 
00650 //--------------------------------------------------------------------------
00651 // void plmap(void (*mapform)(PLINT, PLFLT *, PLFLT *), const char *type,
00652 //            PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat);
00653 //
00654 // plot continental outline in world coordinates
00655 //
00656 // v1.4: machine independant version
00657 // v1.3: replaced plcontinent by plmap, added plmeridians
00658 // v1.2: 2 arguments:  mapform, type of plot
00659 //
00660 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
00661 // coordinate longitudes and latitudes to a plot coordinate system.  By
00662 // using this transform, we can change from a longitude, latitude
00663 // coordinate to a polar stereographic project, for example.  Initially,
00664 // x[0]..[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
00665 // latitudes.  After the call to mapform(), x[] and y[] should be replaced
00666 // by the corresponding plot coordinates.  If no transform is desired,
00667 // mapform can be replaced by NULL.
00668 //
00669 // type is a character string. The value of this parameter determines the
00670 // type of background. The possible values are,
00671 //
00672 //      "globe"         continental outlines
00673 //      "usa"           USA and state boundaries
00674 //      "cglobe"        continental outlines and countries
00675 //      "usaglobe"      USA, state boundaries and continental outlines
00676 // alternatively the filename of a shapefile can be passed if PLplot has
00677 // been compiled with shapelib. In this case either the base name of the
00678 // file can be passed or the filename including the .shp or .shx suffix.
00679 // Only the .shp and .shx files are used.
00680 //
00681 // minlong, maxlong are the values of the longitude on the left and right
00682 // side of the plot, respectively. The value of minlong must be less than
00683 // the values of maxlong, and the values of maxlong-minlong must be less
00684 // or equal to 360.
00685 //
00686 // minlat, maxlat are the minimum and maximum latitudes to be plotted on
00687 // the background.  One can always use -90.0 and 90.0 as the boundary
00688 // outside the plot window will be automatically eliminated.  However, the
00689 // program will be faster if one can reduce the size of the background
00690 // plotted.
00691 //--------------------------------------------------------------------------
00692 
00693 void
00694 plmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type,
00695        PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
00696 {
00697         drawmap( mapform, type, 0.0, 0.0, SHPT_ARC, 0.0, NULL, minlong, maxlong, 
00698                 minlat, maxlat, NULL, 0 );
00699 }
00700 
00701 //--------------------------------------------------------------------------
00702 // void plmapline( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 
00703 //              const char *type, PLFLT minlong, PLFLT maxlong, PLFLT minlat, 
00704 //              PLFLT maxlat, int* plotentries, int nplotentries);
00705 
00706 //New version of plmap which allows us to specify which items in a shapefile
00707 //we want to use. parameters are as above but with the plotentries being an
00708 //array containing the indices of the elements we wish to draw and 
00709 //nplotentries being the number of items in plotentries.
00710 //If shapefile access was not built into plplot then plotentries and 
00711 //nplotentries are ignored. If plotentries is null than all entries are
00712 //drawn and nplotentries is ignored.
00713 //The name distiguishes it from other functions which plot points/text and
00714 //polygons, but note that the type of element in the shapefile need not 
00715 //match the type of element drawn - i.e. arc elements from a shapefile could
00716 //be drawn as points using the plmaptex function.
00717 //--------------------------------------------------------------------------
00718 void
00719 plmapline( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type,
00720        PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, 
00721            int* plotentries, int nplotentries)
00722 {
00723         drawmap( mapform, type, 0.0, 0.0, SHPT_ARC, 0.0, "", minlong, maxlong, 
00724                 minlat, maxlat, plotentries, nplotentries );
00725 }
00726 
00727 //--------------------------------------------------------------------------
00728 // void plmapstring( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 
00729 //              const char *type, PLFLT just, const char *string,
00730 //              PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, 
00731 //              int* plotentries, int nplotentries);
00732 //
00733 //As per plmapline but plots symbols. The map equivalent of plstring. string 
00734 //has the same meaning as in plstring.
00735 //-------------------------------------------------------------------------
00736 void 
00737 plmapstring( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 
00738                 const char *type, const char *string,
00739                 PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, 
00740                 int* plotentries, int nplotentries)
00741 {
00742         drawmap( mapform, type, 1.0, 0.0, SHPT_POINT, 0.5, string, minlong, maxlong,
00743                 minlat, maxlat, plotentries, nplotentries );
00744 }
00745 
00746 //--------------------------------------------------------------------------
00747 // void plmaptex( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 
00748 //              const char *type, PLFLT dx, PLFLT dy PLFLT just, const char *text,
00749 //              PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, 
00750 //              int* plotentries, int nplotentries);
00751 //
00752 //As per plmapline but plots text. The map equivalent of plptex. dx, dy, 
00753 //just and text have the same meaning as in plptex.
00754 //-------------------------------------------------------------------------
00755 void 
00756 plmaptex( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 
00757                 const char *type, PLFLT dx, PLFLT dy, PLFLT just, const char *text,
00758                 PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, 
00759                 int plotentry)
00760 {
00761         drawmap( mapform, type, dx, dy, SHPT_POINT, just, text, minlong, maxlong,
00762                 minlat, maxlat, &plotentry, 1 );
00763 }
00764 
00765 //--------------------------------------------------------------------------
00766 // void plmapfill( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 
00767 //              const char *type, PLFLT minlong, PLFLT maxlong, PLFLT minlat, 
00768 //              PLFLT maxlat, int* plotentries, int nplotentries);
00769 //
00770 //As per plmapline but plots a filled polygon. The map equivalent to 
00771 //plfill. Uses the pattern defined by plsty or plpat.
00772 //-------------------------------------------------------------------------
00773 void 
00774 plmapfill( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 
00775                 const char *type, PLFLT minlong, PLFLT maxlong, PLFLT minlat, 
00776                 PLFLT maxlat, int* plotentries, int nplotentries)
00777 {
00778         drawmap( mapform, type, 0.0, 0.0, SHPT_POLYGON, 0.0, NULL, minlong, maxlong,
00779                 minlat, maxlat, plotentries, nplotentries );
00780 }
00781 
00782 //--------------------------------------------------------------------------
00783 // void plmeridians(void (*mapform)(PLINT, PLFLT *, PLFLT *),
00784 //                  PLFLT dlong, PLFLT dlat, PLFLT minlong, PLFLT maxlong,
00785 //                  PLFLT minlat, PLFLT maxlat);
00786 //
00787 // Plot the latitudes and longitudes on the background.  The lines
00788 // are plotted in the current color and line style.
00789 //
00790 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
00791 // coordinate longitudes and latitudes to a plot coordinate system.  By
00792 // using this transform, we can change from a longitude, latitude
00793 // coordinate to a polar stereographic project, for example.  Initially,
00794 // x[0]..x[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
00795 // latitudes.  After the call to mapform(), x[] and y[] should be replaced
00796 // by the corresponding plot coordinates.  If no transform is desired,
00797 // mapform can be replaced by NULL.
00798 //
00799 // dlat, dlong are the interval in degrees that the latitude and longitude
00800 // lines are to be plotted.
00801 //
00802 // minlong, maxlong are the values of the longitude on the left and right
00803 // side of the plot, respectively. The value of minlong must be less than
00804 // the values of maxlong, and the values of maxlong-minlong must be less
00805 // or equal to 360.
00806 //
00807 // minlat, maxlat are the minimum and maximum latitudes to be plotted on
00808 // the background.  One can always use -90.0 and 90.0 as the boundary
00809 // outside the plot window will be automatically eliminated.  However, the
00810 // program will be faster if one can reduce the size of the background
00811 // plotted.
00812 //--------------------------------------------------------------------------
00813 
00814 #define NSEG    100
00815 
00816 void
00817 plmeridians( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
00818              PLFLT dlong, PLFLT dlat,
00819              PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
00820 {
00821     PLFLT yy, xx, temp, x[2], y[2], dx, dy;
00822 
00823     if ( minlong > maxlong )
00824     {
00825         temp    = minlong;
00826         minlong = maxlong;
00827         maxlong = temp;
00828     }
00829     if ( minlat > maxlat )
00830     {
00831         temp   = minlat;
00832         minlat = maxlat;
00833         maxlat = temp;
00834     }
00835     dx = ( maxlong - minlong ) / NSEG;
00836     dy = ( maxlat - minlat ) / NSEG;
00837 
00838     // latitudes
00839 
00840     for ( yy = dlat * ceil( minlat / dlat ); yy <= maxlat; yy += dlat )
00841     {
00842         if ( mapform == NULL )
00843         {
00844             plpath( NSEG, minlong, yy, maxlong, yy );
00845         }
00846         else
00847         {
00848             for ( xx = minlong; xx < maxlong; xx += dx )
00849             {
00850                 y[0] = y[1] = yy;
00851                 x[0] = xx;
00852                 x[1] = xx + dx;
00853                 ( *mapform )( 2, x, y );
00854                 plline( 2, x, y );
00855             }
00856         }
00857     }
00858 
00859     // longitudes
00860 
00861     for ( xx = dlong * ceil( minlong / dlong ); xx <= maxlong; xx += dlong )
00862     {
00863         if ( mapform == NULL )
00864         {
00865             plpath( NSEG, xx, minlat, xx, maxlat );
00866         }
00867         else
00868         {
00869             for ( yy = minlat; yy < maxlat; yy += dy )
00870             {
00871                 x[0] = x[1] = xx;
00872                 y[0] = yy;
00873                 y[1] = yy + dy;
00874                 ( *mapform )( 2, x, y );
00875                 plline( 2, x, y );
00876             }
00877         }
00878     }
00879 }
00880 
00881 //--------------------------------------------------------------------------
00882 // SHPHandle OpenShapeFile(fn)
00883 //
00897 //--------------------------------------------------------------------------
00898 #ifdef HAVE_SHAPELIB
00899 #ifdef HAVE_SAHOOKS
00900 // Our thanks to Frank Warmerdam, the developer of shapelib for suggesting
00901 // this approach for quieting shapelib "Unable to open" error messages.
00902 static
00903 void CustomErrors( const char *message )
00904 {
00905     if ( strstr( message, "Unable to open" ) == NULL )
00906         fprintf( stderr, "%s\n", message );
00907 }
00908 #endif
00909 
00910 SHPHandle
00911 OpenShapeFile( const char *fn )
00912 {
00913     SHPHandle file;
00914     char      *fs = NULL, *dn = NULL;
00915 #ifdef HAVE_SAHOOKS
00916     SAHooks   sHooks;
00917 
00918     SASetupDefaultHooks( &sHooks );
00919     sHooks.Error = CustomErrors;
00920 #else
00921     // Using ancient version of shapelib without SAHooks or SHPOpenLL.
00922     // For this case live with the misleading "Unable to open" error
00923     // messages.
00924     // int sHooks;
00925 #define SHPOpenLL( a, b, c )    SHPOpen( a, b )
00926 #endif
00927 
00928 //***   search build tree               ***
00929 
00930     if ( plInBuildTree() == 1 )
00931     {
00932         plGetName( SOURCE_DIR, "data", fn, &fs );
00933 
00934         if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
00935             goto done;
00936     }
00937 
00938 //***   search PLPLOT_LIB_ENV = $(PLPLOT_LIB)   ***
00939 
00940 #if defined ( PLPLOT_LIB_ENV )
00941     if ( ( dn = getenv( PLPLOT_LIB_ENV ) ) != NULL )
00942     {
00943         plGetName( dn, "", fn, &fs );
00944 
00945         if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
00946             goto done;
00947         fprintf( stderr, PLPLOT_LIB_ENV "=\"%s\"\n", dn ); // what IS set?
00948     }
00949 #endif  // PLPLOT_LIB_ENV
00950 
00951 //***   search current directory        ***
00952 
00953     if ( ( file = SHPOpenLL( fn, "rb", &sHooks ) ) != NULL )
00954     {
00955         pldebug( "OpenShapeFile", "Found file %s in current directory.\n", fn );
00956         free_mem( fs );
00957         return ( file );
00958     }
00959 
00960 //***   search PLPLOT_HOME_ENV/lib = $(PLPLOT_HOME)/lib ***
00961 
00962 #if defined ( PLPLOT_HOME_ENV )
00963     if ( ( dn = getenv( PLPLOT_HOME_ENV ) ) != NULL )
00964     {
00965         plGetName( dn, "lib", fn, &fs );
00966 
00967         if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
00968             goto done;
00969         fprintf( stderr, PLPLOT_HOME_ENV "=\"%s\"\n", dn ); // what IS set?
00970     }
00971 #endif  // PLPLOT_HOME_ENV/lib
00972 
00973 //***   search installed location       ***
00974 
00975 #if defined ( DATA_DIR )
00976     plGetName( DATA_DIR, "", fn, &fs );
00977 
00978     if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
00979         goto done;
00980 #endif  // DATA_DIR
00981 
00982 //***   search hardwired location       ***
00983 
00984 #ifdef PLLIBDEV
00985     plGetName( PLLIBDEV, "", fn, &fs );
00986 
00987     if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
00988         goto done;
00989 #endif  // PLLIBDEV
00990 
00991 //***   not found, give up      ***
00992     pldebug( "OpenShapeFile", "File %s not found.\n", fn );
00993     free_mem( fs );
00994     return NULL;
00995 
00996 done:
00997     pldebug( "OpenShapeFile", "SHPOpen successfully opened two files with basename %s\n", fs );
00998     free_mem( fs );
00999     return ( file );
01000 }
01001 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines