PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
plmap.c
Go to the documentation of this file.
1 // Continental Outline and Political Boundary Backgrounds
2 //
3 // Some plots need a geographical background such as the global
4 // surface temperatures or the population density. The routine
5 // plmap() will draw one of the following backgrounds: continental
6 // outlines, political boundaries, the United States, and the United
7 // States with the continental outlines. The routine plmeridians()
8 // will add the latitudes and longitudes to the background. After
9 // the background has been drawn, one can use a contour routine or a
10 // symbol plotter to finish off the plot.
11 //
12 // Copyright (C) 1991, 1993, 1994 Wesley Ebisuzaki
13 // Copyright (C) 1994, 2000, 2001 Maurice LeBrun
14 // Copyright (C) 1999 Geoffrey Furnish
15 // Copyright (C) 2000-2014 Alan W. Irwin
16 // Copyright (C) 2001 Andrew Roach
17 // Copyright (C) 2001, 2004 Rafael Laboissiere
18 // Copyright (C) 2002 Vincent Darley
19 // Copyright (C) 2003 Joao Cardoso
20 //
21 // This file is part of PLplot.
22 //
23 // PLplot is free software; you can redistribute it and/or modify
24 // it under the terms of the GNU Library General Public License
25 // as published by the Free Software Foundation; version 2 of the
26 // License.
27 //
28 // PLplot is distributed in the hope that it will be useful, but
29 // WITHOUT ANY WARRANTY; without even the implied warranty of
30 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 // GNU Library General Public License for more details.
32 //
33 // You should have received a copy of the GNU Library General Public
34 // License along with this library; if not, write to the Free Software
35 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301
36 // USA
37 //
38 //
39 
40 #define DEBUG
41 #define NEED_PLDEBUG
42 
43 #include "plplotP.h"
44 
45 #ifdef HAVE_SHAPELIB
46 #include <shapefil.h>
47 
48 SHPHandle
49 OpenShapeFile( const char *fn );
50 
51 #ifdef HAVE_SAHOOKS
52 static void
53 CustomErrors( const char *message );
54 #endif //HAVE_SAHOOKS
55 
56 #define MAP_FILE ""
57 #define OpenMap OpenShapeFile
58 #define CloseMap SHPClose
59 
60 #else //HAVE_SHAPELIB
61 
62 #define MAP_FILE ".map"
63 #define OpenMap plLibOpenPdfstrm
64 #define CloseMap pdf_close
65 #define OFFSET ( 180 * 100 )
66 #define SCALE 100.0
67 #define W_BUFSIZ ( 32 * 1024 )
68 #define SHPT_ARC 1
69 #define SHPT_POINT 2
70 #define SHPT_POLYGON 3
71 
72 #endif //HAVE_SHAPELIB
73 
74 //redistributes the lon value onto either 0-360 or -180-180 for wrapping
75 //purposes.
76 void
77 rebaselon( PLFLT *lon, PLFLT midlon )
78 {
79  if( *lon > midlon + 180.0 )
80  *lon -= floor( ( *lon - midlon - 180.0 ) / 360.0 + 1.0 ) * 360.0;
81  else if( *lon < midlon - 180.0 )
82  *lon += floor( ( midlon - 180.0 - *lon ) / 360.0 + 1.0 ) * 360.0;
83 
84 }
85 
86 //append a PLFLT to an array of PLFLTs. array is a pointer to the array,
87 //n is the current size of the array, val is the value to append.
88 //returns 0 for success, 1 for could not allocate new memory. If memory
89 //could not be allocated then the previous array remains intact
90 int appendflt( PLFLT **array, size_t n, PLFLT val )
91 {
92  size_t i;
93  PLFLT *temp = ( PLFLT* ) malloc( ( n + 1 ) * sizeof( PLFLT ) );
94  if( !temp)
95  return 1;
96  for( i = 0; i < n; ++i )
97  temp[i]=(*array)[i];
98  temp[n] = val;
99  free ( *array );
100  *array = temp;
101  return 0;
102 }
103 
104 //As for appendflt, but for an array of ints
105 int appendint( int **array, size_t n, int val )
106 {
107  size_t i;
108  int *temp = ( int* ) malloc( ( n + 1 ) * sizeof( int ) );
109  if( !temp)
110  return 1;
111  for( i = 0; i < n; ++i )
112  temp[i]=(*array)[i];
113  temp[n] = val;
114  free ( *array );
115  *array = temp;
116  return 0;
117 }
118 
119 //As for appendflt, but for an array of pointers to PLFLTs
120 int appendfltptr( PLFLT ***array, size_t n, PLFLT *val )
121 {
122  size_t i;
123  PLFLT **temp = ( PLFLT** ) malloc( ( n + 1 ) * sizeof( PLFLT *) );
124  if( !temp)
125  return 1;
126  for(i = 0; i < n; ++i)
127  temp[i]=(*array)[i];
128  temp[n] = val;
129  free ( *array );
130  *array = temp;
131  return 0;
132 }
133 
134 //Check to see if the mapform wraps round longitudes. For example, a polar
135 // projection wraps round longitudes, but a cylindrical projection does not.
136 //Returns 1 if the mapform wraps or 0 if not.
137 char
138 checkwrap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), PLFLT lon, PLFLT lat )
139 {
140  PLFLT x[]={lon};
141  PLFLT y[]={lat};
142  PLFLT resultx;
143  PLFLT resulty;
144 
145  if( !mapform )
146  return 0;
147 
148  mapform( 1, x, y );
149  resultx = x[0];
150  resulty = y[0];
151  x[0] = lon + 360;
152  y[0] = lat;
153  return ( ( ABS( x[0] - resultx ) < 1.0e-12 ) && ( ABS( y[0] - resulty ) < 1.0e-12 ) );
154 }
155 
156 //-------------------------------------------------------------------------------------
157 //Actually draw the map lines points and text.
158 //-------------------------------------------------------------------------------------
159 void
160 drawmapdata( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), int shapetype, PLINT n, PLFLT *x, PLFLT *y, PLFLT dx, PLFLT dy, PLFLT just, const char *text )
161 {
162  PLINT i;
163 
164  //do the transform if needed
165  if ( mapform != NULL )
166  ( *mapform )( n, x, y );
167 
168  if ( shapetype == SHPT_ARC )
169  plline( n, x, y );
170  else if ( shapetype == SHPT_POINT )
171  for ( i = 0; i < n; ++i )
172  plptex( x[i], y[i], dx, dy, just, text );
173  else if ( shapetype == SHPT_POLYGON )
174  plfill( n, x, y );
175 #ifdef HAVE_SHAPELIB
176  else if( shapetype == SHPT_ARCZ || shapetype == SHPT_ARCM )
177  plline( n, x, y );
178  else if( shapetype == SHPT_POLYGON || shapetype == SHPT_POLYGONZ || shapetype == SHPT_POLYGONM )
179  plfill( n, x, y );
180  else if( shapetype == SHPT_POINT || shapetype ==SHPT_POINTM || shapetype == SHPT_POINTZ )
181  for ( i = 0; i < n; ++i )
182  plptex( x[i], y[i], dx, dy, just, text );
183 #endif //HAVE_SHAPELIB
184 }
185 
186 
187 //-------------------------------------------------------------------------------------
188 //This is a function called by the front end map functions to do the map drawing. Its
189 //parameters are:
190 //mapform: The transform used to convert the data in raw coordinates to x, y positions
191 //on the plot
192 //type: either one of the plplot provided lat/lon maps or the path/file name of a
193 //shapefile
194 //dx/dy: the gradient of text/symbols drawn if text is non-null
195 //shapetype: one of ARC, SHPT_ARCZ, SHPT_ARCM, SHPT_POLYGON, SHPT_POLYGONZ,
196 //SHPT_POLYGONM, SHPT_POINT, SHPT_POINTM, SHPT_POINTZ. See drawmapdata() for the
197 //how each type is rendered. But Basically the ARC options are lines, the POLYGON
198 //options are filled polygons, the POINT options are points/text. Options beginning
199 //SHPT will only be defined if HAVE_SHAPELIB is true
200 //text: The text (which can be actual text or a unicode symbol) to be drawn at
201 //each point
202 //minlong/maxlong: The min/max longitude when using a plplot provided map or x value if
203 //using a shapefile
204 //minlat/maxlat: The min/max latitude when using a plplot provided map or y value if
205 //using a shapefile
206 //plotentries: used only for shapefiles, as one shapefile contains multiple vectors
207 //each representing a different item (e.g. multiple boundaries, multiple height
208 //contours etc. plotentries is an array containing the indices of the
209 //entries within the shapefile that you wish to plot. if plotentries is null all
210 //entries are plotted
211 //nplotentries: the number of elements in plotentries. Ignored if plplot was not built
212 //with shapefile support or if plotentries is null
213 //-------------------------------------------------------------------------------------
214 void
215 drawmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type,
216  PLFLT dx, PLFLT dy, int shapetype, PLFLT just, const char *text,
217  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, int* plotentries, int nplotentries )
218 {
219 #if defined ( HAVE_SHAPELIB ) || defined ( PL_DEPRECATED )
220  int i, j;
221  char *filename = NULL;
222  char truncatedfilename[900];
223  char warning[1024];
224  int nVertices = 200;
225  PLFLT minsectlon, maxsectlon, minsectlat, maxsectlat;
226  PLFLT *bufx = NULL, *bufy = NULL;
227  int bufsize = 0;
228  int filenamelen;
229  PLFLT **splitx = NULL;
230  PLFLT **splity = NULL;
231  int *splitsectionlengths = NULL;
232  int nsplitsections;
233  PLFLT lastsplitpointx;
234  PLFLT lastsplitpointy;
235  PLFLT penultimatesplitpointx;
236  PLFLT penultimatesplitpointy;
237  char islatlon = 1;
238 
239 
240 #ifdef HAVE_SHAPELIB
241  SHPHandle in;
242  int nentries;
243  int entryindex=0;
244  // Unnecessarily set nparts to quiet -O3 -Wuninitialized warnings.
245  //int nparts = 0;
246  int entrynumber = 0;
247  int partnumber = 0;
248  double mins[4];
249  double maxs[4];
250  SHPObject *object = NULL;
251  double *bufxraw;
252  double *bufyraw;
253  char *prjfilename = NULL;
254  PDFstrm *prjfile;
255  char prjtype[]={0, 0, 0, 0, 0, 0, 0};
256  int appendresult=0;
257 #else
258  PDFstrm *in;
259  //PLFLT bufx[ncopies][200], bufy[ncopies][200];
260  unsigned char n_buff[2], buff[800];
261  long int t;
262 #endif
263 
264  //
265  // read map outline
266  //
267 
268  //strip the .shp extension if a shapefile has been provided and add
269  //the needed map file extension if we are not using shapefile
270  if( strstr( type, ".shp" ) )
271  filenamelen=( int )( type-strstr( type, ".shp" ) );
272  else
273  filenamelen = ( int )strlen( type );
274  filename = ( char* )malloc( filenamelen + strlen( MAP_FILE ) + 1 );
275  if(!filename)
276  {
277  plabort("Could not allocate memory for concatenating map filename");
278  return;
279  }
280  strncpy( filename, type, filenamelen );
281  filename[ filenamelen ] = '\0';
282  strcat( filename, MAP_FILE );
283 
284  //copy the filename to a fixed length array in case it is needed for warning messages
285  if( strlen( filename ) < 899 )
286  strcpy( truncatedfilename, filename );
287  else
288  {
289  memcpy( truncatedfilename, filename, 896 );
290  truncatedfilename[896] = '.';
291  truncatedfilename[897] = '.';
292  truncatedfilename[898] = '.';
293  truncatedfilename[899] = '\0';
294  }
295 
296  strcpy( warning, "Could not find " );
297  strcat( warning, filename );
298  strcat( warning, " file." );
299 #ifdef HAVE_SHAPELIB
300  //Open the shp and shx file using shapelib
301  if ( ( in = OpenShapeFile( filename ) ) == NULL )
302  {
303  plabort( warning );
304  free(filename);
305  return;
306  }
307  SHPGetInfo( in, &nentries, &shapetype, mins, maxs );
308  //also check for a prj file which will tell us if the data is lat/lon or projected
309  //if it is projected then set ncopies to 1 - i.e. don't wrap round longitudes
310  prjfilename = ( char* )malloc( filenamelen + 5 );
311  if(!prjfilename)
312  {
313  plabort("Could not allocate memory for generating map projection filename");
314  free(filename);
315  return;
316  }
317  strncpy( prjfilename, type, filenamelen );
318  prjfilename[ filenamelen ] = '\0';
319  strcat( prjfilename, ".prj" );
320  prjfile = plLibOpenPdfstrm( prjfilename );
321  if ( prjfile && prjfile->file )
322  {
323  fread( prjtype, 1, 6, prjfile->file );
324  if( strcmp ( prjtype, "PROJCS" ) == 0 )
325  islatlon = 0;
326  pdf_close( prjfile );
327  }
328  free ( prjfilename );
329  prjfilename = NULL;
330 #else
331  if ( ( in = plLibOpenPdfstrm( filename ) ) == NULL )
332  {
333  plwarn( warning );
334  return;
335  }
336 #endif
337 
338  bufx = NULL;
339  bufy = NULL;
340 
341  for (;; )
342  {
343 #ifdef HAVE_SHAPELIB
344  //each object in the shapefile is split into parts.
345  //If we are need to plot the first part of an object then read in a new object
346  //and check how many parts it has. Otherwise use the object->panPartStart vector
347  //to check the offset of this part and the next part and allocate memory. Copy
348  //the data to this memory converting it to PLFLT and draw it.
349  //finally increment the part number or if we have finished with the object reset the
350  //part numberand increment the object.
351 
352  //break condition if we've reached the end of the file
353  if ( ( !plotentries && ( entrynumber == nentries ) ) || ( plotentries && ( entryindex == nplotentries ) ) )
354  break;
355 
356  //if partnumber == 0 then we need to load the next object
357  if ( partnumber == 0 )
358  {
359  if( plotentries )
360  object = SHPReadObject( in, plotentries[entryindex] );
361  else
362  object = SHPReadObject( in, entrynumber );
363  }
364  //if the object could not be read, increment the object index to read and
365  //return to the top of the loop to try the next object.
366  if ( object == NULL )
367  {
368  entrynumber++;
369  entryindex++;
370  partnumber = 0;
371  continue;
372  }
373 
374  //work out how many points are in the current part
375  if( object->nParts == 0 )
376  nVertices = object->nVertices; //if object->nParts==0, we can still have 1 vertex. A bit odd but it's the way it goes
377  else if ( partnumber == ( object->nParts - 1 ) )
378  nVertices = object->nVertices - object->panPartStart[partnumber];//panPartStart holds the offset for each part
379  else
380  nVertices = object->panPartStart[partnumber + 1] - object->panPartStart[partnumber];//panPartStart holds the offset for each part
381 #endif
382  //allocate memory for the data
383  if ( nVertices > bufsize )
384  {
385  bufsize = nVertices;
386  free( bufx );
387  free( bufy );
388  bufx = (PLFLT*)malloc( (size_t) bufsize * sizeof ( PLFLT ) );
389  bufy = (PLFLT*)malloc( (size_t) bufsize * sizeof ( PLFLT ) );
390  if(!bufx || ! bufy)
391  {
392  plabort("Could not allocate memory for map data");
393  free(filename);
394  free (bufx);
395  free (bufy);
396  return;
397  }
398  }
399 
400 #ifdef HAVE_SHAPELIB
401  //point the plot buffer to the correct starting vertex
402  //and copy it to the PLFLT arrays. If we had object->nParts == 0
403  //then panPartStart will be NULL
404  if(object->nParts>0)
405  {
406  bufxraw = object->padfX + object->panPartStart[partnumber];
407  bufyraw = object->padfY + object->panPartStart[partnumber];
408  }
409  else
410  {
411  bufxraw = object->padfX;
412  bufyraw = object->padfY;
413  }
414 
415  for ( i = 0; i < nVertices; i++ )
416  {
417  bufx[i] = (PLFLT) bufxraw[i];
418  bufy[i] = (PLFLT) bufyraw[i];
419  }
420 
421  //set the min x/y of the object
422  minsectlon = object->dfXMin;
423  maxsectlon = object->dfXMax;
424  minsectlat = object->dfYMin;
425  maxsectlat = object->dfYMax;
426 
427  //increment the partnumber or if we've reached the end of
428  //an entry increment the entrynumber and set partnumber to 0
429  if ( partnumber == object->nParts - 1 || object->nParts == 0 )
430  {
431  entrynumber++;
432  entryindex++;
433  partnumber = 0;
434  SHPDestroyObject( object );
435  object=NULL;
436  }
437  else
438  partnumber++;
439 
440  if (nVertices==0)
441  continue;
442 #else
443  // read in # points in segment
444  if ( pdf_rdx( n_buff, (long) sizeof ( unsigned char ) * 2, in ) == 0 )
445  break;
446  nVertices = ( n_buff[0] << 8 ) + n_buff[1];
447  if ( nVertices == 0 )
448  break;
449 
450  pdf_rdx( buff, (long) sizeof ( unsigned char ) * 4 * nVertices, in );
451  if ( nVertices == 1 )
452  continue;
453 
454  for ( j = i = 0; i < nVertices; i++, j += 2 )
455  {
456  t = ( buff[j] << 8 ) + buff[j + 1];
457  bufx[i] = ( (PLFLT) t - OFFSET ) / SCALE;
458  }
459  for ( i = 0; i < nVertices; i++, j += 2 )
460  {
461  t = ( buff[j] << 8 ) + buff[j + 1];
462  bufy[i] = ( (PLFLT) t - OFFSET ) / SCALE;
463  }
464  //set the min/max section lat/lon with extreme values
465  //to be overwritten later
466  minsectlon = 1000.;
467  maxsectlon = -1000.;
468  minsectlat = 1000.;
469  maxsectlat = -1000.;
470 
471 #endif
472 
473  if( islatlon )
474  {
475  //two obvious issues exist here with plotting longitudes:
476  //
477  //1) wraparound causing lines which go the wrong way round
478  // the globe
479  //2) some people plot lon from 0-360 deg, others from -180 - +180
480  //
481  //we can cure these problems by conditionally adding/subtracting
482  //360 degrees to each data point in order to ensure that the
483  //distance between adgacent points is always less than 180
484  //degrees, then plotting up to 2 out of 5 copies of the data
485  //each separated by 360 degrees.
486 
487  //arrays of pointers to the starts of each section of data that
488  //has been split due to longitude wrapping, and an array of ints
489  //to hold their lengths. Start with splitx and splity having one
490  //element pointing to the beginning of bufx and bufy
491  splitx = ( PLFLT** ) malloc ( sizeof ( PLFLT* ) );
492  splity = ( PLFLT** ) malloc ( sizeof ( PLFLT* ) );
493  //lengths of the split sections
494  splitsectionlengths = ( int* ) malloc ( sizeof ( size_t ) );
495  if(!splitx || ! splity || ! splitsectionlengths)
496  {
497  plabort("Could not allocate memory for longitudinally split map data");
498  free( filename );
499  free ( bufx );
500  free ( bufy );
501  free ( splitx );
502  free ( splity );
503  free ( splitsectionlengths );
504  return;
505  }
506  splitsectionlengths[0]=nVertices;
507  nsplitsections=1;
508  splitx[0]=bufx;
509  splity[0]=bufy;
510 
511  //set the min/max lats/lons
512  minsectlon = MIN( minsectlon, bufx[0] );
513  maxsectlon = MAX( minsectlon, bufx[0] );
514  minsectlat = MIN( minsectlat, bufy[0] );
515  maxsectlat = MAX( maxsectlat, bufy[0] );
516  //ensure our lat and lon are on 0-360 grid and split the
517  //data where it wraps.
518  rebaselon( &bufx[0], ( minlong + maxlong ) / 2.0 );
519  for ( i = 1; i < nVertices; i++ )
520  {
521  //put lon into 0-360 degree range
522  rebaselon( &bufx[i], ( minlong + maxlong ) / 2.0 );
523 
524  //check if the previous point is more than 180 degrees away
525  if ( bufx[i - 1] - bufx[i] > 180. || bufx[i - 1] - bufx[i] < -180.)
526  {
527  //check if the map transform deals with wrapping itself, e.g. in a polar projection
528  //in this case give one point overlap to the sections so that lines are contiguous
529  if( checkwrap( mapform, bufx[i], bufy[i] ) )
530  {
531  appendresult += appendfltptr( &splitx, nsplitsections, bufx+i );
532  appendresult += appendfltptr( &splity, nsplitsections, bufy+i );
533  appendresult += appendint( &splitsectionlengths, nsplitsections, nVertices - i );
534  splitsectionlengths[nsplitsections-1] -= splitsectionlengths[nsplitsections] - 1;
535  nsplitsections++;
536  }
537  //if the transform doesn't deal with wrapping then allow 2 points overlap to fill in the
538  //edges
539  else
540  {
541  appendresult += appendfltptr( &splitx, nsplitsections, bufx+i-1 );
542  appendresult += appendfltptr( &splity, nsplitsections, bufy+i-1 );
543  appendresult += appendint( &splitsectionlengths, nsplitsections, nVertices - i + 1 );
544  splitsectionlengths[nsplitsections-1] -= splitsectionlengths[nsplitsections] - 2;
545  nsplitsections++;
546  }
547  if( appendresult > 0 )
548  {
549  plabort("Could not allocate memory for appending to longitudinally split map data");
550  free( filename );
551  free ( bufx );
552  free ( bufy );
553  free ( splitx );
554  free ( splity );
555  free ( splitsectionlengths );
556  return;
557  }
558 
559  }
560 
561  //update the mins and maxs
562  minsectlon = MIN( minsectlon, bufx[i] );
563  maxsectlon = MAX( minsectlon, bufx[i] );
564  minsectlat = MIN( minsectlat, bufy[i] );
565  maxsectlat = MAX( maxsectlat, bufy[i] );
566  }
567 
568  //check if the latitude and longitude range means we need to plot this section
569  if ( ( maxsectlat > minlat ) && ( minsectlat < maxlat )
570  && ( maxsectlon > minlong ) && ( minsectlon < maxlong ))
571  {
572  //plot each split in turn, now is where we deal with the end points to
573  //ensure we draw to the edge of the map
574  for( i = 0; i < nsplitsections; ++i )
575  {
576  //check if the first 2 or last 1 points of the split section need
577  //wrapping and add or subtract 360 from them. Note that when the next
578  //section is drawn the code below will undo this if needed
579  if( splitsectionlengths[i] > 2 )
580  {
581  if ( splitx[i][1] - splitx[i][2] > 180. )
582  splitx[i][1] -= 360.0;
583  else if ( splitx[i][1] - splitx[i][2] < -180. )
584  splitx[i][1] += 360.0;
585  }
586 
587  if ( splitx[i][0] - splitx[i][1] > 180. )
588  splitx[i][0] -= 360.0;
589  else if ( splitx[i][0] - splitx[i][1] < -180. )
590  splitx[i][0] += 360.0;
591 
592  if ( splitx[i][splitsectionlengths[i] - 2] - splitx[i][splitsectionlengths[i] - 1] > 180. )
593  splitx[i][splitsectionlengths[i] - 1] += 360.0;
594  else if ( splitx[i][splitsectionlengths[i] - 2] - splitx[i][splitsectionlengths[i] - 1] < -180. )
595  splitx[i][splitsectionlengths[i] - 1] -= 360.0;
596 
597  //save the last 2 points - they will be needed by the next
598  //split section and will be overwritten by the mapform
599  lastsplitpointx = splitx[i][splitsectionlengths[i]-1];
600  lastsplitpointy = splity[i][splitsectionlengths[i]-1];
601  penultimatesplitpointx = splitx[i][splitsectionlengths[i]-2];
602  penultimatesplitpointy = splity[i][splitsectionlengths[i]-2];
603 
604  //draw the split section
605  drawmapdata(mapform, shapetype, splitsectionlengths[i], splitx[i], splity[i], dx, dy, just, text );
606 
607  for(j=1; j<splitsectionlengths[i]; ++j)
608  {
609  if((splitx[i][j]<200.0 && splitx[i][j-1] > 260.0) || ( splitx[i][j-1]<200.0 && splitx[i][j] > 260.0 ))
610  plwarn("wrapping error");
611  }
612 
613  //restore the last 2 points
614  splitx[i][splitsectionlengths[i]-1] = lastsplitpointx;
615  splity[i][splitsectionlengths[i]-1] = lastsplitpointy;
616  splitx[i][splitsectionlengths[i]-2] = penultimatesplitpointx;
617  splity[i][splitsectionlengths[i]-2] = penultimatesplitpointy;
618  }
619  }
620  }
621  else
622  {
623  drawmapdata(mapform, shapetype, nVertices, bufx, bufy, dx, dy, just, text );
624  }
625 
626  free( splitx );
627  free( splity );
628  free( splitsectionlengths );
629 
630 
631  }
632  // Close map file
633 #ifdef HAVE_SHAPELIB
634  SHPClose( in );
635 #else
636  pdf_close( in );
637 #endif
638 
639  //free memory
640  free( bufx );
641  free( bufy );
642  free( filename );
643 #else // defined (HAVE_SHAPELIB) || defined (PL_DEPRECATED)
644  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" );
645 #endif // defined (HAVE_SHAPELIB) || defined (PL_DEPRECATED)
646 }
647 
648 
649 
650 //--------------------------------------------------------------------------
651 // void plmap(void (*mapform)(PLINT, PLFLT *, PLFLT *), const char *type,
652 // PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat);
653 //
654 // plot continental outline in world coordinates
655 //
656 // v1.4: machine independant version
657 // v1.3: replaced plcontinent by plmap, added plmeridians
658 // v1.2: 2 arguments: mapform, type of plot
659 //
660 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
661 // coordinate longitudes and latitudes to a plot coordinate system. By
662 // using this transform, we can change from a longitude, latitude
663 // coordinate to a polar stereographic project, for example. Initially,
664 // x[0]..[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
665 // latitudes. After the call to mapform(), x[] and y[] should be replaced
666 // by the corresponding plot coordinates. If no transform is desired,
667 // mapform can be replaced by NULL.
668 //
669 // type is a character string. The value of this parameter determines the
670 // type of background. The possible values are,
671 //
672 // "globe" continental outlines
673 // "usa" USA and state boundaries
674 // "cglobe" continental outlines and countries
675 // "usaglobe" USA, state boundaries and continental outlines
676 // alternatively the filename of a shapefile can be passed if PLplot has
677 // been compiled with shapelib. In this case either the base name of the
678 // file can be passed or the filename including the .shp or .shx suffix.
679 // Only the .shp and .shx files are used.
680 //
681 // minlong, maxlong are the values of the longitude on the left and right
682 // side of the plot, respectively. The value of minlong must be less than
683 // the values of maxlong, and the values of maxlong-minlong must be less
684 // or equal to 360.
685 //
686 // minlat, maxlat are the minimum and maximum latitudes to be plotted on
687 // the background. One can always use -90.0 and 90.0 as the boundary
688 // outside the plot window will be automatically eliminated. However, the
689 // program will be faster if one can reduce the size of the background
690 // plotted.
691 //--------------------------------------------------------------------------
692 
693 void
694 plmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type,
695  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
696 {
697  drawmap( mapform, type, 0.0, 0.0, SHPT_ARC, 0.0, NULL, minlong, maxlong,
698  minlat, maxlat, NULL, 0 );
699 }
700 
701 //--------------------------------------------------------------------------
702 // void plmapline( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
703 // const char *type, PLFLT minlong, PLFLT maxlong, PLFLT minlat,
704 // PLFLT maxlat, int* plotentries, int nplotentries);
705 
706 //New version of plmap which allows us to specify which items in a shapefile
707 //we want to use. parameters are as above but with the plotentries being an
708 //array containing the indices of the elements we wish to draw and
709 //nplotentries being the number of items in plotentries.
710 //If shapefile access was not built into plplot then plotentries and
711 //nplotentries are ignored. If plotentries is null than all entries are
712 //drawn and nplotentries is ignored.
713 //The name distiguishes it from other functions which plot points/text and
714 //polygons, but note that the type of element in the shapefile need not
715 //match the type of element drawn - i.e. arc elements from a shapefile could
716 //be drawn as points using the plmaptex function.
717 //--------------------------------------------------------------------------
718 void
719 plmapline( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type,
720  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat,
721  int* plotentries, int nplotentries)
722 {
723  drawmap( mapform, type, 0.0, 0.0, SHPT_ARC, 0.0, "", minlong, maxlong,
724  minlat, maxlat, plotentries, nplotentries );
725 }
726 
727 //--------------------------------------------------------------------------
728 // void plmapstring( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
729 // const char *type, PLFLT just, const char *string,
730 // PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat,
731 // int* plotentries, int nplotentries);
732 //
733 //As per plmapline but plots symbols. The map equivalent of plstring. string
734 //has the same meaning as in plstring.
735 //-------------------------------------------------------------------------
736 void
737 plmapstring( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
738  const char *type, const char *string,
739  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat,
740  int* plotentries, int nplotentries)
741 {
742  drawmap( mapform, type, 1.0, 0.0, SHPT_POINT, 0.5, string, minlong, maxlong,
743  minlat, maxlat, plotentries, nplotentries );
744 }
745 
746 //--------------------------------------------------------------------------
747 // void plmaptex( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
748 // const char *type, PLFLT dx, PLFLT dy PLFLT just, const char *text,
749 // PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat,
750 // int* plotentries, int nplotentries);
751 //
752 //As per plmapline but plots text. The map equivalent of plptex. dx, dy,
753 //just and text have the same meaning as in plptex.
754 //-------------------------------------------------------------------------
755 void
756 plmaptex( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
757  const char *type, PLFLT dx, PLFLT dy, PLFLT just, const char *text,
758  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat,
759  int plotentry)
760 {
761  drawmap( mapform, type, dx, dy, SHPT_POINT, just, text, minlong, maxlong,
762  minlat, maxlat, &plotentry, 1 );
763 }
764 
765 //--------------------------------------------------------------------------
766 // void plmapfill( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
767 // const char *type, PLFLT minlong, PLFLT maxlong, PLFLT minlat,
768 // PLFLT maxlat, int* plotentries, int nplotentries);
769 //
770 //As per plmapline but plots a filled polygon. The map equivalent to
771 //plfill. Uses the pattern defined by plsty or plpat.
772 //-------------------------------------------------------------------------
773 void
774 plmapfill( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
775  const char *type, PLFLT minlong, PLFLT maxlong, PLFLT minlat,
776  PLFLT maxlat, int* plotentries, int nplotentries)
777 {
778  drawmap( mapform, type, 0.0, 0.0, SHPT_POLYGON, 0.0, NULL, minlong, maxlong,
779  minlat, maxlat, plotentries, nplotentries );
780 }
781 
782 //--------------------------------------------------------------------------
783 // void plmeridians(void (*mapform)(PLINT, PLFLT *, PLFLT *),
784 // PLFLT dlong, PLFLT dlat, PLFLT minlong, PLFLT maxlong,
785 // PLFLT minlat, PLFLT maxlat);
786 //
787 // Plot the latitudes and longitudes on the background. The lines
788 // are plotted in the current color and line style.
789 //
790 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
791 // coordinate longitudes and latitudes to a plot coordinate system. By
792 // using this transform, we can change from a longitude, latitude
793 // coordinate to a polar stereographic project, for example. Initially,
794 // x[0]..x[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
795 // latitudes. After the call to mapform(), x[] and y[] should be replaced
796 // by the corresponding plot coordinates. If no transform is desired,
797 // mapform can be replaced by NULL.
798 //
799 // dlat, dlong are the interval in degrees that the latitude and longitude
800 // lines are to be plotted.
801 //
802 // minlong, maxlong are the values of the longitude on the left and right
803 // side of the plot, respectively. The value of minlong must be less than
804 // the values of maxlong, and the values of maxlong-minlong must be less
805 // or equal to 360.
806 //
807 // minlat, maxlat are the minimum and maximum latitudes to be plotted on
808 // the background. One can always use -90.0 and 90.0 as the boundary
809 // outside the plot window will be automatically eliminated. However, the
810 // program will be faster if one can reduce the size of the background
811 // plotted.
812 //--------------------------------------------------------------------------
813 
814 #define NSEG 100
815 
816 void
817 plmeridians( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
818  PLFLT dlong, PLFLT dlat,
819  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
820 {
821  PLFLT yy, xx, temp, x[2], y[2], dx, dy;
822 
823  if ( minlong > maxlong )
824  {
825  temp = minlong;
826  minlong = maxlong;
827  maxlong = temp;
828  }
829  if ( minlat > maxlat )
830  {
831  temp = minlat;
832  minlat = maxlat;
833  maxlat = temp;
834  }
835  dx = ( maxlong - minlong ) / NSEG;
836  dy = ( maxlat - minlat ) / NSEG;
837 
838  // latitudes
839 
840  for ( yy = dlat * ceil( minlat / dlat ); yy <= maxlat; yy += dlat )
841  {
842  if ( mapform == NULL )
843  {
844  plpath( NSEG, minlong, yy, maxlong, yy );
845  }
846  else
847  {
848  for ( xx = minlong; xx < maxlong; xx += dx )
849  {
850  y[0] = y[1] = yy;
851  x[0] = xx;
852  x[1] = xx + dx;
853  ( *mapform )( 2, x, y );
854  plline( 2, x, y );
855  }
856  }
857  }
858 
859  // longitudes
860 
861  for ( xx = dlong * ceil( minlong / dlong ); xx <= maxlong; xx += dlong )
862  {
863  if ( mapform == NULL )
864  {
865  plpath( NSEG, xx, minlat, xx, maxlat );
866  }
867  else
868  {
869  for ( yy = minlat; yy < maxlat; yy += dy )
870  {
871  x[0] = x[1] = xx;
872  y[0] = yy;
873  y[1] = yy + dy;
874  ( *mapform )( 2, x, y );
875  plline( 2, x, y );
876  }
877  }
878  }
879 }
880 
881 //--------------------------------------------------------------------------
882 // SHPHandle OpenShapeFile(fn)
883 //
897 //--------------------------------------------------------------------------
898 #ifdef HAVE_SHAPELIB
899 #ifdef HAVE_SAHOOKS
900 // Our thanks to Frank Warmerdam, the developer of shapelib for suggesting
901 // this approach for quieting shapelib "Unable to open" error messages.
902 static
903 void CustomErrors( const char *message )
904 {
905  if ( strstr( message, "Unable to open" ) == NULL )
906  fprintf( stderr, "%s\n", message );
907 }
908 #endif
909 
910 SHPHandle
911 OpenShapeFile( const char *fn )
912 {
913  SHPHandle file;
914  char *fs = NULL, *dn = NULL;
915 #ifdef HAVE_SAHOOKS
916  SAHooks sHooks;
917 
918  SASetupDefaultHooks( &sHooks );
919  sHooks.Error = CustomErrors;
920 #else
921  // Using ancient version of shapelib without SAHooks or SHPOpenLL.
922  // For this case live with the misleading "Unable to open" error
923  // messages.
924  // int sHooks;
925 #define SHPOpenLL( a, b, c ) SHPOpen( a, b )
926 #endif
927 
928 //*** search build tree ***
929 
930  if ( plInBuildTree() == 1 )
931  {
932  plGetName( SOURCE_DIR, "data", fn, &fs );
933 
934  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
935  goto done;
936  }
937 
938 //*** search PLPLOT_LIB_ENV = $(PLPLOT_LIB) ***
939 
940 #if defined ( PLPLOT_LIB_ENV )
941  if ( ( dn = getenv( PLPLOT_LIB_ENV ) ) != NULL )
942  {
943  plGetName( dn, "", fn, &fs );
944 
945  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
946  goto done;
947  fprintf( stderr, PLPLOT_LIB_ENV "=\"%s\"\n", dn ); // what IS set?
948  }
949 #endif // PLPLOT_LIB_ENV
950 
951 //*** search current directory ***
952 
953  if ( ( file = SHPOpenLL( fn, "rb", &sHooks ) ) != NULL )
954  {
955  pldebug( "OpenShapeFile", "Found file %s in current directory.\n", fn );
956  free_mem( fs );
957  return ( file );
958  }
959 
960 //*** search PLPLOT_HOME_ENV/lib = $(PLPLOT_HOME)/lib ***
961 
962 #if defined ( PLPLOT_HOME_ENV )
963  if ( ( dn = getenv( PLPLOT_HOME_ENV ) ) != NULL )
964  {
965  plGetName( dn, "lib", fn, &fs );
966 
967  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
968  goto done;
969  fprintf( stderr, PLPLOT_HOME_ENV "=\"%s\"\n", dn ); // what IS set?
970  }
971 #endif // PLPLOT_HOME_ENV/lib
972 
973 //*** search installed location ***
974 
975 #if defined ( DATA_DIR )
976  plGetName( DATA_DIR, "", fn, &fs );
977 
978  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
979  goto done;
980 #endif // DATA_DIR
981 
982 //*** search hardwired location ***
983 
984 #ifdef PLLIBDEV
985  plGetName( PLLIBDEV, "", fn, &fs );
986 
987  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
988  goto done;
989 #endif // PLLIBDEV
990 
991 //*** not found, give up ***
992  pldebug( "OpenShapeFile", "File %s not found.\n", fn );
993  free_mem( fs );
994  return NULL;
995 
996 done:
997  pldebug( "OpenShapeFile", "SHPOpen successfully opened two files with basename %s\n", fs );
998  free_mem( fs );
999  return ( file );
1000 }
1001 #endif