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