PLplot  5.11.1
 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 //name: 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 //minx/maxx: The min/max longitude when using a plplot provided map or x value if
202 //using a shapefile
203 //miny/maxy: 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 *name,
215  PLFLT dx, PLFLT dy, int shapetype, PLFLT just, const char *text,
216  PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy, const PLINT *plotentries, PLINT 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( name, ".shp" ) )
270  filenamelen = (int) ( name - strstr( name, ".shp" ) );
271  else
272  filenamelen = (int) strlen( name );
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, name, 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  free( filename );
313  plabort( "Could not allocate memory for generating map projection filename" );
314  return;
315  }
316  strncpy( prjfilename, name, 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], ( minx + maxx ) / 2.0 );
518  for ( i = 1; i < nVertices; i++ )
519  {
520  //put lon into 0-360 degree range
521  rebaselon( &bufx[i], ( minx + maxx ) / 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 > miny ) && ( minsectlat < maxy )
568  && ( maxsectlon > minx ) && ( minsectlon < maxx ) )
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 *name,
648 // PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy);
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 // name is a character string. The value of this parameter determines the
666 // name 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 // minx, maxx are the minimum and maximum untransformed x values to be
678 // plotted. For the built in plots these are longitudes. For other
679 //shapefiles these are in the units of the shapefile. The value of minx
680 //must be less than the values of maxx.
681 //
682 // miny, maxy are as per minx and maxx except for the built in plots
683 //the units are latitude.
684 //--------------------------------------------------------------------------
685 
686 void
687 plmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *name,
688  PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy )
689 {
690  drawmap( mapform, name, 0.0, 0.0, SHPT_ARC, 0.0, NULL, minx, maxx,
691  miny, maxy, NULL, 0 );
692 }
693 
694 //--------------------------------------------------------------------------
695 // void plmapline( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
696 // const char *name, PLFLT minx, PLFLT maxx, PLFLT miny,
697 // PLFLT maxy, const PLINT *plotentries, PLINT nplotentries);
698 
699 //New version of plmap which allows us to specify which items in a shapefile
700 //we want to use. parameters are as above but with the plotentries being an
701 //array containing the indices of the elements we wish to draw and
702 //nplotentries being the number of items in plotentries.
703 //If shapefile access was not built into plplot then plotentries and
704 //nplotentries are ignored. If plotentries is null than all entries are
705 //drawn and nplotentries is ignored.
706 //The name distiguishes it from other functions which plot points/text and
707 //polygons, but note that the type of element in the shapefile need not
708 //match the type of element drawn - i.e. arc elements from a shapefile could
709 //be drawn as points using the plmaptex function.
710 //--------------------------------------------------------------------------
711 void
712 plmapline( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *name,
713  PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy,
714  const PLINT *plotentries, PLINT nplotentries )
715 {
716  drawmap( mapform, name, 0.0, 0.0, SHPT_ARC, 0.0, "", minx, maxx,
717  miny, maxy, plotentries, nplotentries );
718 }
719 
720 //--------------------------------------------------------------------------
721 // void plmapstring( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
722 // const char *name, PLFLT just, const char *string,
723 // PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy,
724 // const PLINT *plotentries, PLINT nplotentries);
725 //
726 //As per plmapline but plots symbols. The map equivalent of plstring. string
727 //has the same meaning as in plstring.
728 //--------------------------------------------------------------------------
729 void
730 plmapstring( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
731  const char *name, const char *string,
732  PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy,
733  const PLINT *plotentries, PLINT nplotentries )
734 {
735  drawmap( mapform, name, 1.0, 0.0, SHPT_POINT, 0.5, string, minx, maxx,
736  miny, maxy, plotentries, nplotentries );
737 }
738 
739 //--------------------------------------------------------------------------
740 // void plmaptex( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
741 // const char *name, PLFLT dx, PLFLT dy PLFLT just, const char *text,
742 // PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy,
743 // PLINT plotentry);
744 //
745 //As per plmapline but plots text. The map equivalent of plptex. dx, dy,
746 //just and text have the same meaning as in plptex.
747 //--------------------------------------------------------------------------
748 void
749 plmaptex( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
750  const char *name, PLFLT dx, PLFLT dy, PLFLT just, const char *text,
751  PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy,
752  PLINT plotentry )
753 {
754  drawmap( mapform, name, dx, dy, SHPT_POINT, just, text, minx, maxx,
755  miny, maxy, &plotentry, 1 );
756 }
757 
758 //--------------------------------------------------------------------------
759 // void plmapfill( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
760 // const char *name, PLFLT minx, PLFLT maxx, PLFLT miny,
761 // PLFLT maxy, const PLINT *plotentries, PLINT nplotentries);
762 //
763 //As per plmapline but plots a filled polygon. The map equivalent to
764 //plfill. Uses the pattern defined by plsty or plpat.
765 //--------------------------------------------------------------------------
766 void
767 plmapfill( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
768  const char *name, PLFLT minx, PLFLT maxx, PLFLT miny,
769  PLFLT maxy, const PLINT *plotentries, PLINT nplotentries )
770 {
771  drawmap( mapform, name, 0.0, 0.0, SHPT_POLYGON, 0.0, NULL, minx, maxx,
772  miny, maxy, plotentries, nplotentries );
773 }
774 
775 //--------------------------------------------------------------------------
776 // void plmeridians(void (*mapform)(PLINT, PLFLT *, PLFLT *),
777 // PLFLT dlong, PLFLT dlat, PLFLT minx, PLFLT maxx,
778 // PLFLT miny, PLFLT maxy);
779 //
780 // Plot the latitudes and longitudes on the background. The lines
781 // are plotted in the current color and line style.
782 //
783 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
784 // coordinate longitudes and latitudes to a plot coordinate system. By
785 // using this transform, we can change from a longitude, latitude
786 // coordinate to a polar stereographic project, for example. Initially,
787 // x[0]..x[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
788 // latitudes. After the call to mapform(), x[] and y[] should be replaced
789 // by the corresponding plot coordinates. If no transform is desired,
790 // mapform can be replaced by NULL.
791 //
792 // dlat, dlong are the interval in degrees that the latitude and longitude
793 // lines are to be plotted.
794 //
795 // minx, maxx are the values of the longitude on the left and right
796 // side of the plot, respectively. The value of minx must be less than
797 // the values of maxx, and the values of maxx-minx must be less
798 // or equal to 360.
799 //
800 // miny, maxy are the minimum and maximum latitudes to be plotted on
801 // the background. One can always use -90.0 and 90.0 as the boundary
802 // outside the plot window will be automatically eliminated. However, the
803 // program will be faster if one can reduce the size of the background
804 // plotted.
805 //--------------------------------------------------------------------------
806 
807 #define NSEG 100
808 
809 void
810 plmeridians( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
811  PLFLT dlong, PLFLT dlat,
812  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
813 {
814  PLFLT yy, xx, temp, x[2], y[2], dx, dy;
815 
816  if ( minlong > maxlong )
817  {
818  temp = minlong;
819  minlong = maxlong;
820  maxlong = temp;
821  }
822  if ( minlat > maxlat )
823  {
824  temp = minlat;
825  minlat = maxlat;
826  maxlat = temp;
827  }
828  dx = ( maxlong - minlong ) / NSEG;
829  dy = ( maxlat - minlat ) / NSEG;
830 
831  // latitudes
832 
833  for ( yy = dlat * ceil( minlat / dlat ); yy <= maxlat; yy += dlat )
834  {
835  if ( mapform == NULL )
836  {
837  plpath( NSEG, minlong, yy, maxlong, yy );
838  }
839  else
840  {
841  for ( xx = minlong; xx < maxlong; xx += dx )
842  {
843  y[0] = y[1] = yy;
844  x[0] = xx;
845  x[1] = xx + dx;
846  ( *mapform )( 2, x, y );
847  plline( 2, x, y );
848  }
849  }
850  }
851 
852  // longitudes
853 
854  for ( xx = dlong * ceil( minlong / dlong ); xx <= maxlong; xx += dlong )
855  {
856  if ( mapform == NULL )
857  {
858  plpath( NSEG, xx, minlat, xx, maxlat );
859  }
860  else
861  {
862  for ( yy = minlat; yy < maxlat; yy += dy )
863  {
864  x[0] = x[1] = xx;
865  y[0] = yy;
866  y[1] = yy + dy;
867  ( *mapform )( 2, x, y );
868  plline( 2, x, y );
869  }
870  }
871  }
872 }
873 
874 //--------------------------------------------------------------------------
875 // SHPHandle OpenShapeFile(fn)
876 //
890 //--------------------------------------------------------------------------
891 #ifdef HAVE_SHAPELIB
892 #ifdef HAVE_SAHOOKS
893 // Our thanks to Frank Warmerdam, the developer of shapelib for suggesting
894 // this approach for quieting shapelib "Unable to open" error messages.
895 static
896 void CustomErrors( const char *message )
897 {
898  if ( strstr( message, "Unable to open" ) == NULL )
899  fprintf( stderr, "%s\n", message );
900 }
901 #endif
902 
903 SHPHandle
904 OpenShapeFile( const char *fn )
905 {
906  SHPHandle file;
907  char *fs = NULL, *dn = NULL;
908 #ifdef HAVE_SAHOOKS
909  SAHooks sHooks;
910 
911  SASetupDefaultHooks( &sHooks );
912  sHooks.Error = CustomErrors;
913 #else
914  // Using ancient version of shapelib without SAHooks or SHPOpenLL.
915  // For this case live with the misleading "Unable to open" error
916  // messages.
917  // int sHooks;
918 #define SHPOpenLL( a, b, c ) SHPOpen( a, b )
919 #endif
920 
921 //*** search build tree ***
922 
923  if ( plInBuildTree() == 1 )
924  {
925  plGetName( SOURCE_DIR, "data", fn, &fs );
926 
927  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
928  goto done;
929  }
930 
931 //*** search PLPLOT_LIB_ENV = $(PLPLOT_LIB) ***
932 
933 #if defined ( PLPLOT_LIB_ENV )
934  if ( ( dn = getenv( PLPLOT_LIB_ENV ) ) != NULL )
935  {
936  plGetName( dn, "", fn, &fs );
937 
938  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
939  goto done;
940  fprintf( stderr, PLPLOT_LIB_ENV "=\"%s\"\n", dn ); // what IS set?
941  }
942 #endif // PLPLOT_LIB_ENV
943 
944 //*** search current directory ***
945 
946  if ( ( file = SHPOpenLL( fn, "rb", &sHooks ) ) != NULL )
947  {
948  pldebug( "OpenShapeFile", "Found file %s in current directory.\n", fn );
949  free_mem( fs );
950  return ( file );
951  }
952 
953 //*** search PLPLOT_HOME_ENV/lib = $(PLPLOT_HOME)/lib ***
954 
955 #if defined ( PLPLOT_HOME_ENV )
956  if ( ( dn = getenv( PLPLOT_HOME_ENV ) ) != NULL )
957  {
958  plGetName( dn, "lib", fn, &fs );
959 
960  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
961  goto done;
962  fprintf( stderr, PLPLOT_HOME_ENV "=\"%s\"\n", dn ); // what IS set?
963  }
964 #endif // PLPLOT_HOME_ENV/lib
965 
966 //*** search installed location ***
967 
968 #if defined ( DATA_DIR )
969  plGetName( DATA_DIR, "", fn, &fs );
970 
971  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
972  goto done;
973 #endif // DATA_DIR
974 
975 //*** search hardwired location ***
976 
977 #ifdef PLLIBDEV
978  plGetName( PLLIBDEV, "", fn, &fs );
979 
980  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
981  goto done;
982 #endif // PLLIBDEV
983 
984 //*** not found, give up ***
985  pldebug( "OpenShapeFile", "File %s not found.\n", fn );
986  free_mem( fs );
987  return NULL;
988 
989 done:
990  pldebug( "OpenShapeFile", "SHPOpen successfully opened two files with basename %s\n", fs );
991  free_mem( fs );
992  return ( file );
993 }
994 #endif