PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
plmap.c
Go to the documentation of this file.
1 // $Id: plmap.c 12827 2013-12-09 13:20:01Z andrewross $
2 //
3 // Continental Outline and Political Boundary Backgrounds
4 //
5 // Some plots need a geographical background such as the global
6 // surface temperatures or the population density. The routine
7 // plmap() will draw one of the following backgrounds: continental
8 // outlines, political boundaries, the United States, and the United
9 // States with the continental outlines. The routine plmeridians()
10 // will add the latitudes and longitudes to the background. After
11 // the background has been drawn, one can use a contour routine or a
12 // symbol plotter to finish off the plot.
13 //
14 // Copyright (C) 1991, 1993, 1994 Wesley Ebisuzaki
15 // Copyright (C) 1994, 2000, 2001 Maurice LeBrun
16 // Copyright (C) 1999 Geoffrey Furnish
17 // Copyright (C) 2000, 2001, 2002 Alan W. Irwin
18 // Copyright (C) 2001 Andrew Roach
19 // Copyright (C) 2001, 2004 Rafael Laboissiere
20 // Copyright (C) 2002 Vincent Darley
21 // Copyright (C) 2003 Joao Cardoso
22 //
23 // This file is part of PLplot.
24 //
25 // PLplot is free software; you can redistribute it and/or modify
26 // it under the terms of the GNU Library General Public License
27 // as published by the Free Software Foundation; version 2 of the
28 // License.
29 //
30 // PLplot is distributed in the hope that it will be useful, but
31 // WITHOUT ANY WARRANTY; without even the implied warranty of
32 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33 // GNU Library General Public License for more details.
34 //
35 // You should have received a copy of the GNU Library General Public
36 // License along with this library; if not, write to the Free Software
37 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301
38 // USA
39 //
40 //
41 
42 #define DEBUG
43 #define NEED_PLDEBUG
44 
45 #include "plplotP.h"
46 
47 #ifdef HAVE_SHAPELIB
48 #include <shapefil.h>
49 
50 SHPHandle
51 OpenShapeFile( const char *fn );
52 
53 #ifdef HAVE_SAHOOKS
54 static void
55 CustomErrors( const char *message );
56 #endif
57 
58 #endif
59 
60 //--------------------------------------------------------------------------
61 // void plmap(void (*mapform)(PLINT, PLFLT *, PLFLT *), const char *type,
62 // PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat);
63 //
64 // plot continental outline in world coordinates
65 //
66 // v1.4: machine independant version
67 // v1.3: replaced plcontinent by plmap, added plmeridians
68 // v1.2: 2 arguments: mapform, type of plot
69 //
70 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
71 // coordinate longitudes and latitudes to a plot coordinate system. By
72 // using this transform, we can change from a longitude, latitude
73 // coordinate to a polar stereographic project, for example. Initially,
74 // x[0]..[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
75 // latitudes. After the call to mapform(), x[] and y[] should be replaced
76 // by the corresponding plot coordinates. If no transform is desired,
77 // mapform can be replaced by NULL.
78 //
79 // type is a character string. The value of this parameter determines the
80 // type of background. The possible values are,
81 //
82 // "globe" continental outlines
83 // "usa" USA and state boundaries
84 // "cglobe" continental outlines and countries
85 // "usaglobe" USA, state boundaries and continental outlines
86 // alternatively the filename of a shapefile can be passed if PLplot has
87 // been compiled with shapelib. In this case either the base name of the
88 // file can be passed or the filename including the .shp or .shx suffix.
89 // Only the .shp and .shx files are used.
90 //
91 // minlong, maxlong are the values of the longitude on the left and right
92 // side of the plot, respectively. The value of minlong must be less than
93 // the values of maxlong, and the values of maxlong-minlong must be less
94 // or equal to 360.
95 //
96 // minlat, maxlat are the minimum and maximum latitudes to be plotted on
97 // the background. One can always use -90.0 and 90.0 as the boundary
98 // outside the plot window will be automatically eliminated. However, the
99 // program will be faster if one can reduce the size of the background
100 // plotted.
101 //--------------------------------------------------------------------------
102 
103 #ifdef HAVE_SHAPELIB
104 #define MAP_FILE ""
105 #define OpenMap OpenShapeFile
106 #define CloseMap SHPClose
107 #else
108 #define MAP_FILE ".map"
109 #define OpenMap plLibOpenPdfstrm
110 #define CloseMap pdf_close
111 #define OFFSET ( 180 * 100 )
112 #define SCALE 100.0
113 #define W_BUFSIZ ( 32 * 1024 )
114 #endif
115 
116 void
117 plmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type,
118  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
119 {
120 #if defined ( HAVE_SHAPELIB ) || defined ( PL_DEPRECATED )
121  int i, j;
122  char *filename;
123  char *warning;
124  int n = 200;
125  PLFLT minsectlon, maxsectlon, minsectlat, maxsectlat;
126  const int ncopies = 5; //must be odd - original plus copies shifted by multiples of +/- 360
127  const int mid = ncopies / 2 + 1;
128  PLFLT **bufx = NULL, **bufy = NULL;
129  int bufsize = 0;
130 
131 #ifdef HAVE_SHAPELIB
132  SHPHandle in;
133  int nentries;
134  // Unnecessarily set nparts to quiet -O3 -Wuninitialized warnings.
135  int nparts = 0;
136  int entrynumber = 0;
137  int partnumber = 0;
138  int shapetype;
139  double mins[4];
140  double maxs[4];
141  SHPObject *object = NULL;
142  double *bufxraw;
143  double *bufyraw;
144 #else
145  register PDFstrm *in;
146  //PLFLT bufx[ncopies][200], bufy[ncopies][200];
147  unsigned char n_buff[2], buff[800];
148  long int t;
149  int k;
150 #endif
151 
152  //
153  // read map outline
154  //
155  filename = malloc( strlen( type ) + strlen( MAP_FILE ) + 1 );
156  strcpy( filename, type );
157  strcat( filename, MAP_FILE );
158 
159  warning = malloc( strlen( type ) + strlen( MAP_FILE ) + 50 );
160  strcpy( warning, "Could not find " );
161  strcat( warning, filename );
162  strcat( warning, " file." );
163 #ifdef HAVE_SHAPELIB
164  if ( ( in = OpenShapeFile( filename ) ) == NULL )
165  {
166  plwarn( warning );
167  return;
168  }
169  SHPGetInfo( in, &nentries, &shapetype, mins, maxs );
170 #else
171  if ( ( in = plLibOpenPdfstrm( filename ) ) == NULL )
172  {
173  plwarn( warning );
174  return;
175  }
176 #endif
177 
178  bufx = malloc( (size_t) ncopies * sizeof ( PLFLT* ) );
179  bufy = malloc( (size_t) ncopies * sizeof ( PLFLT* ) );
180  for ( i = 0; i < ncopies; i++ )
181  {
182  bufx[i] = NULL;
183  bufy[i] = NULL;
184  }
185 
186  for (;; )
187  {
188 #ifdef HAVE_SHAPELIB
189  //break condition if we've reached the end of the file
190  if ( entrynumber == nentries )
191  break;
192  //if partnumber == 0 then we need to load the next object
193  if ( partnumber == 0 )
194  {
195  object = SHPReadObject( in, entrynumber );
196  nparts = object->nParts;
197  }
198 
199  //work out how many points are in the current part
200  if ( partnumber == ( nparts - 1 ) )
201  n = object->nVertices - object->panPartStart[partnumber];
202  else
203  n = object->panPartStart[partnumber + 1] - object->panPartStart[partnumber];
204 #endif
205  //allocate memory for the data
206  if ( n > bufsize )
207  {
208  bufsize = n;
209  for ( i = 0; i < ncopies; i++ )
210  {
211  if ( bufx[i] )
212  free( bufx[i] );
213  if ( bufy[i] )
214  free( bufy[i] );
215  bufx[i] = malloc( (size_t) bufsize * sizeof ( double ) );
216  bufy[i] = malloc( (size_t) bufsize * sizeof ( double ) );
217  }
218  }
219 
220 #ifdef HAVE_SHAPELIB
221  //point the plot buffer to the correct starting vertex
222  //and copy it to the PLFLT arrays
223  bufxraw = object->padfX + object->panPartStart[partnumber];
224  bufyraw = object->padfY + object->panPartStart[partnumber];
225  for ( i = 0; i < n; i++ )
226  {
227  bufx[mid][i] = (PLFLT) bufxraw[i];
228  for ( j = 0; j < ncopies; j++ )
229  bufy[j][i] = (PLFLT) bufyraw[i];
230  }
231 
232  //set the minlat/lon of the object
233  minsectlon = object->dfXMin;
234  maxsectlon = object->dfXMax;
235  minsectlat = object->dfYMin;
236  maxsectlat = object->dfYMax;
237 
238  //increment the partnumber or if we've reached the end of
239  //an entry increment the entrynumber and set partnumber to 0
240  if ( partnumber == nparts - 1 )
241  {
242  entrynumber++;
243  partnumber = 0;
244  }
245  else
246  partnumber++;
247 #else
248  // read in # points in segment
249  if ( pdf_rdx( n_buff, (long) sizeof ( unsigned char ) * 2, in ) == 0 )
250  break;
251  n = ( n_buff[0] << 8 ) + n_buff[1];
252  if ( n == 0 )
253  break;
254 
255  pdf_rdx( buff, (long) sizeof ( unsigned char ) * 4 * n, in );
256  if ( n == 1 )
257  continue;
258 
259  for ( j = i = 0; i < n; i++, j += 2 )
260  {
261  t = ( buff[j] << 8 ) + buff[j + 1];
262  bufx[mid][i] = ( (PLFLT) t - OFFSET ) / SCALE;
263  }
264  for ( i = 0; i < n; i++, j += 2 )
265  {
266  t = ( buff[j] << 8 ) + buff[j + 1];
267  bufy[0][i] = ( (PLFLT) t - OFFSET ) / SCALE;
268  for ( k = 1; k < ncopies; k++ )
269  bufy[k][i] = bufy[0][i];
270  }
271  //set the min/max section lat/lon with extreme values
272  //to be overwritten later
273  minsectlon = 1000.;
274  maxsectlon = -1000.;
275  minsectlat = 1000.;
276  maxsectlat = -1000.;
277 
278 #endif
279  //two obvious issues exist here with plotting longitudes:
280  //
281  //1) wraparound causing lines which go the wrong way round
282  // the globe
283  //2) some people plot lon from 0-360 deg, others from -180 - +180
284  //
285  //we can cure these problems by conditionally adding/subtracting
286  //360 degrees to each data point in order to ensure that the
287  //distance between adgacent points is always less than 180
288  //degrees, then plotting up to 2 out of 5 copies of the data
289  //each separated by 360 degrees.
290 
291  for ( i = 0; i < n - 1; i++ )
292  {
293  if ( bufx[mid][i] - bufx[mid][i + 1] > 180. )
294  bufx[mid][i + 1] += 360.;
295  else if ( bufx[mid][i] - bufx[mid][i + 1] < -180. )
296  bufx[mid][i + 1] -= 360.;
297  }
298  for ( i = 0; i < n; i++ )
299  {
300  for ( j = 0; j < mid; j++ )
301  bufx[j][i] = bufx[mid][i] + 360. * (PLFLT) ( j - mid );
302  for ( j = mid + 1; j < ncopies; j++ )
303  bufx[j][i] = bufx[mid][i] + 360. * (PLFLT) ( j - mid );
304 #ifndef HAVE_SHAPELIB
305  minsectlon = MIN( minsectlon, bufx[mid][i] );
306  maxsectlon = MAX( minsectlon, bufx[mid][i] );
307  minsectlat = MIN( minsectlat, bufy[mid][i] );
308  maxsectlat = MAX( minsectlat, bufy[mid][i] );
309 #endif
310  }
311 
312  //check if the latitude range means we need to plot this section
313  if ( ( maxsectlat > minlat ) && ( minsectlat < maxlat ) )
314  {
315  //check which of the translated maps fall within the
316  //range and transform and plot them - note more than one
317  //map may be needed due to wrapping
318  for ( j = 0; j < ncopies; j++ )
319  {
320  if ( ( minsectlon + 360. * (PLFLT) ( j - mid ) < maxlong )
321  && ( maxsectlon + 360. * (PLFLT) ( j - mid ) > minlong ) )
322  {
323  if ( mapform != NULL )
324  ( *mapform )( n, bufx[j], bufy[j] );
325  plline( n, bufx[j], bufy[j] );
326  }
327  }
328  }
329 
330 
331 
332 #ifdef HAVE_SHAPELIB
333  if ( partnumber == 0 )
334  SHPDestroyObject( object );
335 #endif
336  }
337  // Close map file
338 #ifdef HAVE_SHAPELIB
339  SHPClose( in );
340 #else
341  pdf_close( in );
342 #endif
343 
344  //free memory
345  for ( i = 0; i < ncopies; i++ )
346  {
347  if ( bufx[i] )
348  free( bufx[i] );
349  if ( bufy[i] )
350  free( bufy[i] );
351  }
352  free( bufx );
353  free( bufy );
354  free( filename );
355  free( warning );
356 #else // defined (HAVE_SHAPELIB) || defined (PL_DEPRECATED)
357  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" );
358 #endif // defined (HAVE_SHAPELIB) || defined (PL_DEPRECATED)
359 }
360 
361 //--------------------------------------------------------------------------
362 // void plmeridians(void (*mapform)(PLINT, PLFLT *, PLFLT *),
363 // PLFLT dlong, PLFLT dlat, PLFLT minlong, PLFLT maxlong,
364 // PLFLT minlat, PLFLT maxlat);
365 //
366 // Plot the latitudes and longitudes on the background. The lines
367 // are plotted in the current color and line style.
368 //
369 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
370 // coordinate longitudes and latitudes to a plot coordinate system. By
371 // using this transform, we can change from a longitude, latitude
372 // coordinate to a polar stereographic project, for example. Initially,
373 // x[0]..x[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
374 // latitudes. After the call to mapform(), x[] and y[] should be replaced
375 // by the corresponding plot coordinates. If no transform is desired,
376 // mapform can be replaced by NULL.
377 //
378 // dlat, dlong are the interval in degrees that the latitude and longitude
379 // lines are to be plotted.
380 //
381 // minlong, maxlong are the values of the longitude on the left and right
382 // side of the plot, respectively. The value of minlong must be less than
383 // the values of maxlong, and the values of maxlong-minlong must be less
384 // or equal to 360.
385 //
386 // minlat, maxlat are the minimum and maximum latitudes to be plotted on
387 // the background. One can always use -90.0 and 90.0 as the boundary
388 // outside the plot window will be automatically eliminated. However, the
389 // program will be faster if one can reduce the size of the background
390 // plotted.
391 //--------------------------------------------------------------------------
392 
393 #define NSEG 100
394 
395 void
396 plmeridians( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
397  PLFLT dlong, PLFLT dlat,
398  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
399 {
400  PLFLT yy, xx, temp, x[2], y[2], dx, dy;
401 
402  if ( minlong > maxlong )
403  {
404  temp = minlong;
405  minlong = maxlong;
406  maxlong = temp;
407  }
408  if ( minlat > maxlat )
409  {
410  temp = minlat;
411  minlat = maxlat;
412  maxlat = temp;
413  }
414  dx = ( maxlong - minlong ) / NSEG;
415  dy = ( maxlat - minlat ) / NSEG;
416 
417  // latitudes
418 
419  for ( yy = dlat * ceil( minlat / dlat ); yy <= maxlat; yy += dlat )
420  {
421  if ( mapform == NULL )
422  {
423  plpath( NSEG, minlong, yy, maxlong, yy );
424  }
425  else
426  {
427  for ( xx = minlong; xx < maxlong; xx += dx )
428  {
429  y[0] = y[1] = yy;
430  x[0] = xx;
431  x[1] = xx + dx;
432  ( *mapform )( 2, x, y );
433  plline( 2, x, y );
434  }
435  }
436  }
437 
438  // longitudes
439 
440  for ( xx = dlong * ceil( minlong / dlong ); xx <= maxlong; xx += dlong )
441  {
442  if ( mapform == NULL )
443  {
444  plpath( NSEG, xx, minlat, xx, maxlat );
445  }
446  else
447  {
448  for ( yy = minlat; yy < maxlat; yy += dy )
449  {
450  x[0] = x[1] = xx;
451  y[0] = yy;
452  y[1] = yy + dy;
453  ( *mapform )( 2, x, y );
454  plline( 2, x, y );
455  }
456  }
457  }
458 }
459 
460 //--------------------------------------------------------------------------
461 // SHPHandle OpenShapeFile(fn)
462 //
476 //--------------------------------------------------------------------------
477 #ifdef HAVE_SHAPELIB
478 #ifdef HAVE_SAHOOKS
479 // Our thanks to Frank Warmerdam, the developer of shapelib for suggesting
480 // this approach for quieting shapelib "Unable to open" error messages.
481 static
482 void CustomErrors( const char *message )
483 {
484  if ( strstr( message, "Unable to open" ) == NULL )
485  fprintf( stderr, "%s\n", message );
486 }
487 #endif
488 
489 SHPHandle
490 OpenShapeFile( const char *fn )
491 {
492  SHPHandle file;
493  char *fs = NULL, *dn = NULL;
494 #ifdef HAVE_SAHOOKS
495  SAHooks sHooks;
496 
497  SASetupDefaultHooks( &sHooks );
498  sHooks.Error = CustomErrors;
499 #else
500  // Using ancient version of shapelib without SAHooks or SHPOpenLL.
501  // For this case live with the misleading "Unable to open" error
502  // messages.
503  // int sHooks;
504 #define SHPOpenLL( a, b, c ) SHPOpen( a, b )
505 #endif
506 
507 //*** search build tree ***
508 
509  if ( plInBuildTree() == 1 )
510  {
511  plGetName( SOURCE_DIR, "data", fn, &fs );
512 
513  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
514  goto done;
515  }
516 
517 //*** search PLPLOT_LIB_ENV = $(PLPLOT_LIB) ***
518 
519 #if defined ( PLPLOT_LIB_ENV )
520  if ( ( dn = getenv( PLPLOT_LIB_ENV ) ) != NULL )
521  {
522  plGetName( dn, "", fn, &fs );
523 
524  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
525  goto done;
526  fprintf( stderr, PLPLOT_LIB_ENV "=\"%s\"\n", dn ); // what IS set?
527  }
528 #endif // PLPLOT_LIB_ENV
529 
530 //*** search current directory ***
531 
532  if ( ( file = SHPOpenLL( fn, "rb", &sHooks ) ) != NULL )
533  {
534  pldebug( "OpenShapeFile", "Found file %s in current directory.\n", fn );
535  free_mem( fs );
536  return ( file );
537  }
538 
539 //*** search PLPLOT_HOME_ENV/lib = $(PLPLOT_HOME)/lib ***
540 
541 #if defined ( PLPLOT_HOME_ENV )
542  if ( ( dn = getenv( PLPLOT_HOME_ENV ) ) != NULL )
543  {
544  plGetName( dn, "lib", fn, &fs );
545 
546  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
547  goto done;
548  fprintf( stderr, PLPLOT_HOME_ENV "=\"%s\"\n", dn ); // what IS set?
549  }
550 #endif // PLPLOT_HOME_ENV/lib
551 
552 //*** search installed location ***
553 
554 #if defined ( DATA_DIR )
555  plGetName( DATA_DIR, "", fn, &fs );
556 
557  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
558  goto done;
559 #endif // DATA_DIR
560 
561 //*** search hardwired location ***
562 
563 #ifdef PLLIBDEV
564  plGetName( PLLIBDEV, "", fn, &fs );
565 
566  if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
567  goto done;
568 #endif // PLLIBDEV
569 
570 //*** not found, give up ***
571  pldebug( "OpenShapeFile", "File %s not found.\n", fn );
572  free_mem( fs );
573  return NULL;
574 
575 done:
576  pldebug( "OpenShapeFile", "SHPOpen successfully opened two files with basename %s\n", fs );
577  free_mem( fs );
578  return ( file );
579 }
580 #endif