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