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