PLplot  5.9.9
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
plmap.c
Go to the documentation of this file.
1 // $Id: plmap.c 12281 2012-11-23 21:27:51Z 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 NEED_PLDEBUG 1
43 
44 #include "plplotP.h"
45 
46 #ifdef HAVE_SHAPELIB
47 #include <shapefil.h>
48 
49 SHPHandle
50 OpenShapeFile( const char *fn );
51 
52 #endif
53 
54 //--------------------------------------------------------------------------
55 // void plmap(void (*mapform)(PLINT, PLFLT *, PLFLT *), const char *type,
56 // PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat);
57 //
58 // plot continental outline in world coordinates
59 //
60 // v1.4: machine independant version
61 // v1.3: replaced plcontinent by plmap, added plmeridians
62 // v1.2: 2 arguments: mapform, type of plot
63 //
64 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
65 // coordinate longitudes and latitudes to a plot coordinate system. By
66 // using this transform, we can change from a longitude, latitude
67 // coordinate to a polar stereographic project, for example. Initially,
68 // x[0]..[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
69 // latitudes. After the call to mapform(), x[] and y[] should be replaced
70 // by the corresponding plot coordinates. If no transform is desired,
71 // mapform can be replaced by NULL.
72 //
73 // type is a character string. The value of this parameter determines the
74 // type of background. The possible values are,
75 //
76 // "globe" continental outlines
77 // "usa" USA and state boundaries
78 // "cglobe" continental outlines and countries
79 // "usaglobe" USA, state boundaries and continental outlines
80 // alternatively the filename of a shapefile can be passed if PLplot has
81 // been compiled with shapelib. In this case either the base name of the
82 // file can be passed or the filename including the .shp or .shx suffix.
83 // Only the .shp and .shx files are used.
84 //
85 // minlong, maxlong are the values of the longitude on the left and right
86 // side of the plot, respectively. The value of minlong must be less than
87 // the values of maxlong, and the values of maxlong-minlong must be less
88 // or equal to 360.
89 //
90 // minlat, maxlat are the minimum and maximum latitudes to be plotted on
91 // the background. One can always use -90.0 and 90.0 as the boundary
92 // outside the plot window will be automatically eliminated. However, the
93 // program will be faster if one can reduce the size of the background
94 // plotted.
95 //--------------------------------------------------------------------------
96 
97 #ifdef HAVE_SHAPELIB
98 #define MAP_FILE ""
99 #define OpenMap OpenShapeFile
100 #define CloseMap SHPClose
101 #else
102 #define MAP_FILE ".map"
103 #define OpenMap plLibOpenPdfstrm
104 #define CloseMap pdf_close
105 #define OFFSET ( 180 * 100 )
106 #define SCALE 100.0
107 #define W_BUFSIZ ( 32 * 1024 )
108 #endif
109 
110 void
111 plmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type,
112  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
113 {
114 #if defined ( HAVE_SHAPELIB ) || defined ( PL_DEPRECATED )
115  int i, j;
116  char *filename;
117  char *warning;
118  int n = 200;
119  PLFLT minsectlon, maxsectlon, minsectlat, maxsectlat;
120  const int ncopies = 5; //must be odd - original plus copies shifted by multiples of +/- 360
121  const int mid = ncopies / 2 + 1;
122  PLFLT **bufx = NULL, **bufy = NULL;
123  int bufsize = 0;
124 
125 #ifdef HAVE_SHAPELIB
126  SHPHandle in;
127  int nentries;
128  int nparts;
129  int entrynumber = 0;
130  int partnumber = 0;
131  int shapetype;
132  double mins[4];
133  double maxs[4];
134  SHPObject *object = NULL;
135  double *bufxraw;
136  double *bufyraw;
137 #else
138  register PDFstrm *in;
139  //PLFLT bufx[ncopies][200], bufy[ncopies][200];
140  unsigned char n_buff[2], buff[800];
141  long int t;
142  int k;
143 #endif
144 
145  //
146  // read map outline
147  //
148  filename = malloc( strlen( type ) + strlen( MAP_FILE ) + 1 );
149  strcpy( filename, type );
150  strcat( filename, MAP_FILE );
151 
152  warning = malloc( strlen( type ) + strlen( MAP_FILE ) + 50 );
153  strcpy( warning, "Could not find " );
154  strcat( warning, filename );
155  strcat( warning, " file." );
156 #ifdef HAVE_SHAPELIB
157  if ( ( in = OpenShapeFile( filename ) ) == NULL )
158  {
159  plwarn( warning );
160  return;
161  }
162  SHPGetInfo( in, &nentries, &shapetype, mins, maxs );
163 #else
164  if ( ( in = plLibOpenPdfstrm( filename ) ) == NULL )
165  {
166  plwarn( warning );
167  return;
168  }
169 #endif
170 
171  bufx = malloc( ncopies * sizeof ( PLFLT* ) );
172  bufy = malloc( ncopies * sizeof ( PLFLT* ) );
173  for ( i = 0; i < ncopies; i++ )
174  {
175  bufx[i] = NULL;
176  bufy[i] = NULL;
177  }
178 
179  for (;; )
180  {
181 #ifdef HAVE_SHAPELIB
182  //break condition if we've reached the end of the file
183  if ( entrynumber == nentries )
184  break;
185  //if partnumber == 0 then we need to load the next object
186  if ( partnumber == 0 )
187  {
188  object = SHPReadObject( in, entrynumber );
189  nparts = object->nParts;
190  }
191 
192  //work out how many points are in the current part
193  if ( partnumber == ( nparts - 1 ) )
194  n = object->nVertices - object->panPartStart[partnumber];
195  else
196  n = object->panPartStart[partnumber + 1] - object->panPartStart[partnumber];
197 #endif
198  //allocate memory for the data
199  if ( n > bufsize )
200  {
201  bufsize = n;
202  for ( i = 0; i < ncopies; i++ )
203  {
204  if ( bufx[i] )
205  free( bufx[i] );
206  if ( bufy[i] )
207  free( bufy[i] );
208  bufx[i] = malloc( bufsize * sizeof ( double ) );
209  bufy[i] = malloc( bufsize * sizeof ( double ) );
210  }
211  }
212 
213 #ifdef HAVE_SHAPELIB
214  //point the plot buffer to the correct starting vertex
215  //and copy it to the PLFLT arrays
216  bufxraw = object->padfX + object->panPartStart[partnumber];
217  bufyraw = object->padfY + object->panPartStart[partnumber];
218  for ( i = 0; i < n; i++ )
219  {
220  bufx[mid][i] = (PLFLT) bufxraw[i];
221  for ( j = 0; j < ncopies; j++ )
222  bufy[j][i] = (PLFLT) bufyraw[i];
223  }
224 
225  //set the minlat/lon of the object
226  minsectlon = object->dfXMin;
227  maxsectlon = object->dfXMax;
228  minsectlat = object->dfYMin;
229  maxsectlat = object->dfYMax;
230 
231  //increment the partnumber or if we've reached the end of
232  //an entry increment the entrynumber and set partnumber to 0
233  if ( partnumber == nparts - 1 )
234  {
235  entrynumber++;
236  partnumber = 0;
237  }
238  else
239  partnumber++;
240 #else
241  // read in # points in segment
242  if ( pdf_rdx( n_buff, (long) sizeof ( unsigned char ) * 2, in ) == 0 )
243  break;
244  n = ( n_buff[0] << 8 ) + n_buff[1];
245  if ( n == 0 )
246  break;
247 
248  pdf_rdx( buff, (long) sizeof ( unsigned char ) * 4 * n, in );
249  if ( n == 1 )
250  continue;
251 
252  for ( j = i = 0; i < n; i++, j += 2 )
253  {
254  t = ( buff[j] << 8 ) + buff[j + 1];
255  bufx[mid][i] = ( (PLFLT) t - OFFSET ) / SCALE;
256  }
257  for ( i = 0; i < n; i++, j += 2 )
258  {
259  t = ( buff[j] << 8 ) + buff[j + 1];
260  bufy[0][i] = ( (PLFLT) t - OFFSET ) / SCALE;
261  for ( k = 1; k < ncopies; k++ )
262  bufy[k][i] = bufy[0][i];
263  }
264  //set the min/max section lat/lon with extreme values
265  //to be overwritten later
266  minsectlon = 1000.;
267  maxsectlon = -1000.;
268  minsectlat = 1000.;
269  maxsectlat = -1000.;
270 
271 #endif
272  //two obvious issues exist here with plotting longitudes:
273  //
274  //1) wraparound causing lines which go the wrong way round
275  // the globe
276  //2) some people plot lon from 0-360 deg, others from -180 - +180
277  //
278  //we can cure these problems by conditionally adding/subtracting
279  //360 degrees to each data point in order to ensure that the
280  //distance between adgacent points is always less than 180
281  //degrees, then plotting up to 2 out of 5 copies of the data
282  //each separated by 360 degrees.
283 
284  for ( i = 0; i < n - 1; i++ )
285  {
286  if ( bufx[mid][i] - bufx[mid][i + 1] > 180. )
287  bufx[mid][i + 1] += 360.;
288  else if ( bufx[mid][i] - bufx[mid][i + 1] < -180. )
289  bufx[mid][i + 1] -= 360.;
290  }
291  for ( i = 0; i < n; i++ )
292  {
293  for ( j = 0; j < mid; j++ )
294  bufx[j][i] = bufx[mid][i] + 360. * (PLFLT) ( j - mid );
295  for ( j = mid + 1; j < ncopies; j++ )
296  bufx[j][i] = bufx[mid][i] + 360. * (PLFLT) ( j - mid );
297 #ifndef HAVE_SHAPELIB
298  minsectlon = MIN( minsectlon, bufx[mid][i] );
299  maxsectlon = MAX( minsectlon, bufx[mid][i] );
300  minsectlat = MIN( minsectlat, bufy[mid][i] );
301  maxsectlat = MAX( minsectlat, bufy[mid][i] );
302 #endif
303  }
304 
305  //check if the latitude range means we need to plot this section
306  if ( ( maxsectlat > minlat ) && ( minsectlat < maxlat ) )
307  {
308  //check which of the translated maps fall within the
309  //range and transform and plot them - note more than one
310  //map may be needed due to wrapping
311  for ( j = 0; j < ncopies; j++ )
312  {
313  if ( ( minsectlon + 360. * (PLFLT) ( j - mid ) < maxlong )
314  && ( maxsectlon + 360. * (PLFLT) ( j - mid ) > minlong ) )
315  {
316  if ( mapform != NULL )
317  ( *mapform )( n, bufx[j], bufy[j] );
318  plline( n, bufx[j], bufy[j] );
319  }
320  }
321  }
322 
323 
324 
325 #ifdef HAVE_SHAPELIB
326  if ( partnumber == 0 )
327  SHPDestroyObject( object );
328 #endif
329  }
330  // Close map file
331 #ifdef HAVE_SHAPELIB
332  SHPClose( in );
333 #else
334  pdf_close( in );
335 #endif
336 
337  //free memory
338  for ( i = 0; i < ncopies; i++ )
339  {
340  if ( bufx[i] )
341  free( bufx[i] );
342  if ( bufy[i] )
343  free( bufy[i] );
344  }
345  free( bufx );
346  free( bufy );
347  free( filename );
348  free( warning );
349 #else // defined (HAVE_SHAPELIB) || defined (PL_DEPRECATED)
350  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" );
351 #endif // defined (HAVE_SHAPELIB) || defined (PL_DEPRECATED)
352 }
353 
354 //--------------------------------------------------------------------------
355 // void plmeridians(void (*mapform)(PLINT, PLFLT *, PLFLT *),
356 // PLFLT dlong, PLFLT dlat, PLFLT minlong, PLFLT maxlong,
357 // PLFLT minlat, PLFLT maxlat);
358 //
359 // Plot the latitudes and longitudes on the background. The lines
360 // are plotted in the current color and line style.
361 //
362 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
363 // coordinate longitudes and latitudes to a plot coordinate system. By
364 // using this transform, we can change from a longitude, latitude
365 // coordinate to a polar stereographic project, for example. Initially,
366 // x[0]..x[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
367 // latitudes. After the call to mapform(), x[] and y[] should be replaced
368 // by the corresponding plot coordinates. If no transform is desired,
369 // mapform can be replaced by NULL.
370 //
371 // dlat, dlong are the interval in degrees that the latitude and longitude
372 // lines are to be plotted.
373 //
374 // minlong, maxlong are the values of the longitude on the left and right
375 // side of the plot, respectively. The value of minlong must be less than
376 // the values of maxlong, and the values of maxlong-minlong must be less
377 // or equal to 360.
378 //
379 // minlat, maxlat are the minimum and maximum latitudes to be plotted on
380 // the background. One can always use -90.0 and 90.0 as the boundary
381 // outside the plot window will be automatically eliminated. However, the
382 // program will be faster if one can reduce the size of the background
383 // plotted.
384 //--------------------------------------------------------------------------
385 
386 #define NSEG 100
387 
388 void
389 plmeridians( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
390  PLFLT dlong, PLFLT dlat,
391  PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
392 {
393  PLFLT yy, xx, temp, x[2], y[2], dx, dy;
394 
395  if ( minlong > maxlong )
396  {
397  temp = minlong;
398  minlong = maxlong;
399  maxlong = temp;
400  }
401  if ( minlat > maxlat )
402  {
403  temp = minlat;
404  minlat = maxlat;
405  maxlat = temp;
406  }
407  dx = ( maxlong - minlong ) / NSEG;
408  dy = ( maxlat - minlat ) / NSEG;
409 
410  // latitudes
411 
412  for ( yy = dlat * ceil( minlat / dlat ); yy <= maxlat; yy += dlat )
413  {
414  if ( mapform == NULL )
415  {
416  plpath( NSEG, minlong, yy, maxlong, yy );
417  }
418  else
419  {
420  for ( xx = minlong; xx < maxlong; xx += dx )
421  {
422  y[0] = y[1] = yy;
423  x[0] = xx;
424  x[1] = xx + dx;
425  ( *mapform )( 2, x, y );
426  plline( 2, x, y );
427  }
428  }
429  }
430 
431  // longitudes
432 
433  for ( xx = dlong * ceil( minlong / dlong ); xx <= maxlong; xx += dlong )
434  {
435  if ( mapform == NULL )
436  {
437  plpath( NSEG, xx, minlat, xx, maxlat );
438  }
439  else
440  {
441  for ( yy = minlat; yy < maxlat; yy += dy )
442  {
443  x[0] = x[1] = xx;
444  y[0] = yy;
445  y[1] = yy + dy;
446  ( *mapform )( 2, x, y );
447  plline( 2, x, y );
448  }
449  }
450  }
451 }
452 
453 //--------------------------------------------------------------------------
454 // SHPHandle OpenShapeFile(fn)
455 //
469 //--------------------------------------------------------------------------
470 #ifdef HAVE_SHAPELIB
471 SHPHandle
472 OpenShapeFile( const char *fn )
473 {
474  SHPHandle file;
475  char *fs = NULL, *dn = NULL;
476 
477 //*** search build tree ***
478 
479  if ( plInBuildTree() == 1 )
480  {
481  plGetName( SOURCE_DIR, "data", fn, &fs );
482 
483  if ( ( file = SHPOpen( fs, "rb" ) ) != NULL )
484  goto done;
485  }
486 
487 //*** search PLPLOT_LIB_ENV = $(PLPLOT_LIB) ***
488 
489 #if defined ( PLPLOT_LIB_ENV )
490  if ( ( dn = getenv( PLPLOT_LIB_ENV ) ) != NULL )
491  {
492  plGetName( dn, "", fn, &fs );
493 
494  if ( ( file = SHPOpen( fs, "rb" ) ) != NULL )
495  goto done;
496  fprintf( stderr, PLPLOT_LIB_ENV "=\"%s\"\n", dn ); // what IS set?
497  }
498 #endif // PLPLOT_LIB_ENV
499 
500 //*** search current directory ***
501 
502  if ( ( file = SHPOpen( fn, "rb" ) ) != NULL )
503  {
504  pldebug( "OpenShapeFile", "Found file %s in current directory.\n", fn );
505  free_mem( fs );
506  return ( file );
507  }
508 
509 //*** search PLPLOT_HOME_ENV/lib = $(PLPLOT_HOME)/lib ***
510 
511 #if defined ( PLPLOT_HOME_ENV )
512  if ( ( dn = getenv( PLPLOT_HOME_ENV ) ) != NULL )
513  {
514  plGetName( dn, "lib", fn, &fs );
515 
516  if ( ( file = SHPOpen( fs, "rb" ) ) != NULL )
517  goto done;
518  fprintf( stderr, PLPLOT_HOME_ENV "=\"%s\"\n", dn ); // what IS set?
519  }
520 #endif // PLPLOT_HOME_ENV/lib
521 
522 //*** search installed location ***
523 
524 #if defined ( DATA_DIR )
525  plGetName( DATA_DIR, "", fn, &fs );
526 
527  if ( ( file = SHPOpen( fs, "rb" ) ) != NULL )
528  goto done;
529 #endif // DATA_DIR
530 
531 //*** search hardwired location ***
532 
533 #ifdef PLLIBDEV
534  plGetName( PLLIBDEV, "", fn, &fs );
535 
536  if ( ( file = SHPOpen( fs, "rb" ) ) != NULL )
537  goto done;
538 #endif // PLLIBDEV
539 
540 //*** not found, give up ***
541  pldebug( "OpenShapeFile", "File %s not found.\n", fn );
542  free_mem( fs );
543  return NULL;
544 
545 done:
546  pldebug( "OpenShapeFile", "Found file %s\n", fs );
547  free_mem( fs );
548  return ( file );
549 }
550 #endif
#define plpath
Definition: plplot.h:656
int pdf_rdx(U_CHAR *x, long nitems, PDFstrm *pdfs)
Definition: pdfutils.c:464
void plGetName(const char *dir, const char *subdir, const char *filename, char **filespec)
Definition: plctrl.c:2435
void mapform(PLINT n, PLFLT *x, PLFLT *y)
Definition: tclAPI.c:3228
Definition: pdf.h:51
#define OFFSET
Definition: plmap.c:105
#define MAX(a, b)
Definition: dsplint.c:28
void PLFLT PLINT PLINT PLFLT x
#define DATA_DIR
Definition: config.h:27
#define MAP_FILE
Definition: plmap.c:102
PDFstrm * plLibOpenPdfstrm(const char *fn)
Definition: plctrl.c:2245
int PLINT
Definition: plplot.h:175
#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:389
void PLFLT PLINT PLINT PLFLT PLFLT y
#define SCALE
Definition: plmap.c:106
#define PLLIBDEV
Definition: plctrl.c:125
void plmap(void(*mapform)(PLINT, PLFLT *, PLFLT *), const char *type, PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat)
Definition: plmap.c:111
#define PLPLOT_HOME_ENV
Definition: plplotP.h:403
#define SOURCE_DIR
Definition: config.h:277
float PLFLT
Definition: plplot.h:159
int pdf_close(PDFstrm *pdfs)
Definition: pdfutils.c:238
void plwarn(const char *errormsg)
Definition: plctrl.c:1846
#define free_mem(a)
Definition: plplotP.h:187
#define NSEG
Definition: plmap.c:386
#define PLPLOT_LIB_ENV
Definition: plplotP.h:401
#define plline
Definition: plplot.h:655
dx
if { $zoomopts($this,1) == 0 } then {
Definition: Plframe.py:613
PLDLLIMPEXP int plInBuildTree()
Definition: plcore.c:2774