PLplot
5.10.0
|
00001 // Continental Outline and Political Boundary Backgrounds 00002 // 00003 // Some plots need a geographical background such as the global 00004 // surface temperatures or the population density. The routine 00005 // plmap() will draw one of the following backgrounds: continental 00006 // outlines, political boundaries, the United States, and the United 00007 // States with the continental outlines. The routine plmeridians() 00008 // will add the latitudes and longitudes to the background. After 00009 // the background has been drawn, one can use a contour routine or a 00010 // symbol plotter to finish off the plot. 00011 // 00012 // Copyright (C) 1991, 1993, 1994 Wesley Ebisuzaki 00013 // Copyright (C) 1994, 2000, 2001 Maurice LeBrun 00014 // Copyright (C) 1999 Geoffrey Furnish 00015 // Copyright (C) 2000-2014 Alan W. Irwin 00016 // Copyright (C) 2001 Andrew Roach 00017 // Copyright (C) 2001, 2004 Rafael Laboissiere 00018 // Copyright (C) 2002 Vincent Darley 00019 // Copyright (C) 2003 Joao Cardoso 00020 // 00021 // This file is part of PLplot. 00022 // 00023 // PLplot is free software; you can redistribute it and/or modify 00024 // it under the terms of the GNU Library General Public License 00025 // as published by the Free Software Foundation; version 2 of the 00026 // License. 00027 // 00028 // PLplot is distributed in the hope that it will be useful, but 00029 // WITHOUT ANY WARRANTY; without even the implied warranty of 00030 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00031 // GNU Library General Public License for more details. 00032 // 00033 // You should have received a copy of the GNU Library General Public 00034 // License along with this library; if not, write to the Free Software 00035 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 00036 // USA 00037 // 00038 // 00039 00040 #define DEBUG 00041 #define NEED_PLDEBUG 00042 00043 #include "plplotP.h" 00044 00045 #ifdef HAVE_SHAPELIB 00046 #include <shapefil.h> 00047 00048 SHPHandle 00049 OpenShapeFile( const char *fn ); 00050 00051 #ifdef HAVE_SAHOOKS 00052 static void 00053 CustomErrors( const char *message ); 00054 #endif //HAVE_SAHOOKS 00055 00056 #define MAP_FILE "" 00057 #define OpenMap OpenShapeFile 00058 #define CloseMap SHPClose 00059 00060 #else //HAVE_SHAPELIB 00061 00062 #define MAP_FILE ".map" 00063 #define OpenMap plLibOpenPdfstrm 00064 #define CloseMap pdf_close 00065 #define OFFSET ( 180 * 100 ) 00066 #define SCALE 100.0 00067 #define W_BUFSIZ ( 32 * 1024 ) 00068 #define SHPT_ARC 1 00069 #define SHPT_POINT 2 00070 #define SHPT_POLYGON 3 00071 00072 #endif //HAVE_SHAPELIB 00073 00074 //redistributes the lon value onto either 0-360 or -180-180 for wrapping 00075 //purposes. 00076 void 00077 rebaselon( PLFLT *lon, PLFLT midlon ) 00078 { 00079 if( *lon > midlon + 180.0 ) 00080 *lon -= floor( ( *lon - midlon - 180.0 ) / 360.0 + 1.0 ) * 360.0; 00081 else if( *lon < midlon - 180.0 ) 00082 *lon += floor( ( midlon - 180.0 - *lon ) / 360.0 + 1.0 ) * 360.0; 00083 00084 } 00085 00086 //append a PLFLT to an array of PLFLTs. array is a pointer to the array, 00087 //n is the current size of the array, val is the value to append. 00088 //returns 0 for success, 1 for could not allocate new memory. If memory 00089 //could not be allocated then the previous array remains intact 00090 int appendflt( PLFLT **array, size_t n, PLFLT val ) 00091 { 00092 size_t i; 00093 PLFLT *temp = ( PLFLT* ) malloc( ( n + 1 ) * sizeof( PLFLT ) ); 00094 if( !temp) 00095 return 1; 00096 for( i = 0; i < n; ++i ) 00097 temp[i]=(*array)[i]; 00098 temp[n] = val; 00099 free ( *array ); 00100 *array = temp; 00101 return 0; 00102 } 00103 00104 //As for appendflt, but for an array of ints 00105 int appendint( int **array, size_t n, int val ) 00106 { 00107 size_t i; 00108 int *temp = ( int* ) malloc( ( n + 1 ) * sizeof( int ) ); 00109 if( !temp) 00110 return 1; 00111 for( i = 0; i < n; ++i ) 00112 temp[i]=(*array)[i]; 00113 temp[n] = val; 00114 free ( *array ); 00115 *array = temp; 00116 return 0; 00117 } 00118 00119 //As for appendflt, but for an array of pointers to PLFLTs 00120 int appendfltptr( PLFLT ***array, size_t n, PLFLT *val ) 00121 { 00122 size_t i; 00123 PLFLT **temp = ( PLFLT** ) malloc( ( n + 1 ) * sizeof( PLFLT *) ); 00124 if( !temp) 00125 return 1; 00126 for(i = 0; i < n; ++i) 00127 temp[i]=(*array)[i]; 00128 temp[n] = val; 00129 free ( *array ); 00130 *array = temp; 00131 return 0; 00132 } 00133 00134 //Check to see if the mapform wraps round longitudes. For example, a polar 00135 // projection wraps round longitudes, but a cylindrical projection does not. 00136 //Returns 1 if the mapform wraps or 0 if not. 00137 char 00138 checkwrap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), PLFLT lon, PLFLT lat ) 00139 { 00140 PLFLT x[]={lon}; 00141 PLFLT y[]={lat}; 00142 PLFLT resultx; 00143 PLFLT resulty; 00144 00145 if( !mapform ) 00146 return 0; 00147 00148 mapform( 1, x, y ); 00149 resultx = x[0]; 00150 resulty = y[0]; 00151 x[0] = lon + 360; 00152 y[0] = lat; 00153 return ( ( ABS( x[0] - resultx ) < 1.0e-12 ) && ( ABS( y[0] - resulty ) < 1.0e-12 ) ); 00154 } 00155 00156 //------------------------------------------------------------------------------------- 00157 //Actually draw the map lines points and text. 00158 //------------------------------------------------------------------------------------- 00159 void 00160 drawmapdata( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), int shapetype, PLINT n, PLFLT *x, PLFLT *y, PLFLT dx, PLFLT dy, PLFLT just, const char *text ) 00161 { 00162 PLINT i; 00163 00164 //do the transform if needed 00165 if ( mapform != NULL ) 00166 ( *mapform )( n, x, y ); 00167 00168 if ( shapetype == SHPT_ARC ) 00169 plline( n, x, y ); 00170 else if ( shapetype == SHPT_POINT ) 00171 for ( i = 0; i < n; ++i ) 00172 plptex( x[i], y[i], dx, dy, just, text ); 00173 else if ( shapetype == SHPT_POLYGON ) 00174 plfill( n, x, y ); 00175 #ifdef HAVE_SHAPELIB 00176 else if( shapetype == SHPT_ARCZ || shapetype == SHPT_ARCM ) 00177 plline( n, x, y ); 00178 else if( shapetype == SHPT_POLYGON || shapetype == SHPT_POLYGONZ || shapetype == SHPT_POLYGONM ) 00179 plfill( n, x, y ); 00180 else if( shapetype == SHPT_POINT || shapetype ==SHPT_POINTM || shapetype == SHPT_POINTZ ) 00181 for ( i = 0; i < n; ++i ) 00182 plptex( x[i], y[i], dx, dy, just, text ); 00183 #endif //HAVE_SHAPELIB 00184 } 00185 00186 00187 //------------------------------------------------------------------------------------- 00188 //This is a function called by the front end map functions to do the map drawing. Its 00189 //parameters are: 00190 //mapform: The transform used to convert the data in raw coordinates to x, y positions 00191 //on the plot 00192 //type: either one of the plplot provided lat/lon maps or the path/file name of a 00193 //shapefile 00194 //dx/dy: the gradient of text/symbols drawn if text is non-null 00195 //shapetype: one of ARC, SHPT_ARCZ, SHPT_ARCM, SHPT_POLYGON, SHPT_POLYGONZ, 00196 //SHPT_POLYGONM, SHPT_POINT, SHPT_POINTM, SHPT_POINTZ. See drawmapdata() for the 00197 //how each type is rendered. But Basically the ARC options are lines, the POLYGON 00198 //options are filled polygons, the POINT options are points/text. Options beginning 00199 //SHPT will only be defined if HAVE_SHAPELIB is true 00200 //text: The text (which can be actual text or a unicode symbol) to be drawn at 00201 //each point 00202 //minlong/maxlong: The min/max longitude when using a plplot provided map or x value if 00203 //using a shapefile 00204 //minlat/maxlat: The min/max latitude when using a plplot provided map or y value if 00205 //using a shapefile 00206 //plotentries: used only for shapefiles, as one shapefile contains multiple vectors 00207 //each representing a different item (e.g. multiple boundaries, multiple height 00208 //contours etc. plotentries is an array containing the indices of the 00209 //entries within the shapefile that you wish to plot. if plotentries is null all 00210 //entries are plotted 00211 //nplotentries: the number of elements in plotentries. Ignored if plplot was not built 00212 //with shapefile support or if plotentries is null 00213 //------------------------------------------------------------------------------------- 00214 void 00215 drawmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type, 00216 PLFLT dx, PLFLT dy, int shapetype, PLFLT just, const char *text, 00217 PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, int* plotentries, int nplotentries ) 00218 { 00219 #if defined ( HAVE_SHAPELIB ) || defined ( PL_DEPRECATED ) 00220 int i, j; 00221 char *filename = NULL; 00222 char truncatedfilename[900]; 00223 char warning[1024]; 00224 int nVertices = 200; 00225 PLFLT minsectlon, maxsectlon, minsectlat, maxsectlat; 00226 PLFLT *bufx = NULL, *bufy = NULL; 00227 int bufsize = 0; 00228 int filenamelen; 00229 PLFLT **splitx = NULL; 00230 PLFLT **splity = NULL; 00231 int *splitsectionlengths = NULL; 00232 int nsplitsections; 00233 PLFLT lastsplitpointx; 00234 PLFLT lastsplitpointy; 00235 PLFLT penultimatesplitpointx; 00236 PLFLT penultimatesplitpointy; 00237 char islatlon = 1; 00238 00239 00240 #ifdef HAVE_SHAPELIB 00241 SHPHandle in; 00242 int nentries; 00243 int entryindex=0; 00244 // Unnecessarily set nparts to quiet -O3 -Wuninitialized warnings. 00245 //int nparts = 0; 00246 int entrynumber = 0; 00247 int partnumber = 0; 00248 double mins[4]; 00249 double maxs[4]; 00250 SHPObject *object = NULL; 00251 double *bufxraw; 00252 double *bufyraw; 00253 char *prjfilename = NULL; 00254 PDFstrm *prjfile; 00255 char prjtype[]={0, 0, 0, 0, 0, 0, 0}; 00256 int appendresult=0; 00257 #else 00258 PDFstrm *in; 00259 //PLFLT bufx[ncopies][200], bufy[ncopies][200]; 00260 unsigned char n_buff[2], buff[800]; 00261 long int t; 00262 #endif 00263 00264 // 00265 // read map outline 00266 // 00267 00268 //strip the .shp extension if a shapefile has been provided and add 00269 //the needed map file extension if we are not using shapefile 00270 if( strstr( type, ".shp" ) ) 00271 filenamelen=( int )( type-strstr( type, ".shp" ) ); 00272 else 00273 filenamelen = ( int )strlen( type ); 00274 filename = ( char* )malloc( filenamelen + strlen( MAP_FILE ) + 1 ); 00275 if(!filename) 00276 { 00277 plabort("Could not allocate memory for concatenating map filename"); 00278 return; 00279 } 00280 strncpy( filename, type, filenamelen ); 00281 filename[ filenamelen ] = '\0'; 00282 strcat( filename, MAP_FILE ); 00283 00284 //copy the filename to a fixed length array in case it is needed for warning messages 00285 if( strlen( filename ) < 899 ) 00286 strcpy( truncatedfilename, filename ); 00287 else 00288 { 00289 memcpy( truncatedfilename, filename, 896 ); 00290 truncatedfilename[896] = '.'; 00291 truncatedfilename[897] = '.'; 00292 truncatedfilename[898] = '.'; 00293 truncatedfilename[899] = '\0'; 00294 } 00295 00296 strcpy( warning, "Could not find " ); 00297 strcat( warning, filename ); 00298 strcat( warning, " file." ); 00299 #ifdef HAVE_SHAPELIB 00300 //Open the shp and shx file using shapelib 00301 if ( ( in = OpenShapeFile( filename ) ) == NULL ) 00302 { 00303 plabort( warning ); 00304 free(filename); 00305 return; 00306 } 00307 SHPGetInfo( in, &nentries, &shapetype, mins, maxs ); 00308 //also check for a prj file which will tell us if the data is lat/lon or projected 00309 //if it is projected then set ncopies to 1 - i.e. don't wrap round longitudes 00310 prjfilename = ( char* )malloc( filenamelen + 5 ); 00311 if(!prjfilename) 00312 { 00313 plabort("Could not allocate memory for generating map projection filename"); 00314 free(filename); 00315 return; 00316 } 00317 strncpy( prjfilename, type, filenamelen ); 00318 prjfilename[ filenamelen ] = '\0'; 00319 strcat( prjfilename, ".prj" ); 00320 prjfile = plLibOpenPdfstrm( prjfilename ); 00321 if ( prjfile && prjfile->file ) 00322 { 00323 fread( prjtype, 1, 6, prjfile->file ); 00324 if( strcmp ( prjtype, "PROJCS" ) == 0 ) 00325 islatlon = 0; 00326 pdf_close( prjfile ); 00327 } 00328 free ( prjfilename ); 00329 prjfilename = NULL; 00330 #else 00331 if ( ( in = plLibOpenPdfstrm( filename ) ) == NULL ) 00332 { 00333 plwarn( warning ); 00334 return; 00335 } 00336 #endif 00337 00338 bufx = NULL; 00339 bufy = NULL; 00340 00341 for (;; ) 00342 { 00343 #ifdef HAVE_SHAPELIB 00344 //each object in the shapefile is split into parts. 00345 //If we are need to plot the first part of an object then read in a new object 00346 //and check how many parts it has. Otherwise use the object->panPartStart vector 00347 //to check the offset of this part and the next part and allocate memory. Copy 00348 //the data to this memory converting it to PLFLT and draw it. 00349 //finally increment the part number or if we have finished with the object reset the 00350 //part numberand increment the object. 00351 00352 //break condition if we've reached the end of the file 00353 if ( ( !plotentries && ( entrynumber == nentries ) ) || ( plotentries && ( entryindex == nplotentries ) ) ) 00354 break; 00355 00356 //if partnumber == 0 then we need to load the next object 00357 if ( partnumber == 0 ) 00358 { 00359 if( plotentries ) 00360 object = SHPReadObject( in, plotentries[entryindex] ); 00361 else 00362 object = SHPReadObject( in, entrynumber ); 00363 } 00364 //if the object could not be read, increment the object index to read and 00365 //return to the top of the loop to try the next object. 00366 if ( object == NULL ) 00367 { 00368 entrynumber++; 00369 entryindex++; 00370 partnumber = 0; 00371 continue; 00372 } 00373 00374 //work out how many points are in the current part 00375 if( object->nParts == 0 ) 00376 nVertices = object->nVertices; //if object->nParts==0, we can still have 1 vertex. A bit odd but it's the way it goes 00377 else if ( partnumber == ( object->nParts - 1 ) ) 00378 nVertices = object->nVertices - object->panPartStart[partnumber];//panPartStart holds the offset for each part 00379 else 00380 nVertices = object->panPartStart[partnumber + 1] - object->panPartStart[partnumber];//panPartStart holds the offset for each part 00381 #endif 00382 //allocate memory for the data 00383 if ( nVertices > bufsize ) 00384 { 00385 bufsize = nVertices; 00386 free( bufx ); 00387 free( bufy ); 00388 bufx = (PLFLT*)malloc( (size_t) bufsize * sizeof ( PLFLT ) ); 00389 bufy = (PLFLT*)malloc( (size_t) bufsize * sizeof ( PLFLT ) ); 00390 if(!bufx || ! bufy) 00391 { 00392 plabort("Could not allocate memory for map data"); 00393 free(filename); 00394 free (bufx); 00395 free (bufy); 00396 return; 00397 } 00398 } 00399 00400 #ifdef HAVE_SHAPELIB 00401 //point the plot buffer to the correct starting vertex 00402 //and copy it to the PLFLT arrays. If we had object->nParts == 0 00403 //then panPartStart will be NULL 00404 if(object->nParts>0) 00405 { 00406 bufxraw = object->padfX + object->panPartStart[partnumber]; 00407 bufyraw = object->padfY + object->panPartStart[partnumber]; 00408 } 00409 else 00410 { 00411 bufxraw = object->padfX; 00412 bufyraw = object->padfY; 00413 } 00414 00415 for ( i = 0; i < nVertices; i++ ) 00416 { 00417 bufx[i] = (PLFLT) bufxraw[i]; 00418 bufy[i] = (PLFLT) bufyraw[i]; 00419 } 00420 00421 //set the min x/y of the object 00422 minsectlon = object->dfXMin; 00423 maxsectlon = object->dfXMax; 00424 minsectlat = object->dfYMin; 00425 maxsectlat = object->dfYMax; 00426 00427 //increment the partnumber or if we've reached the end of 00428 //an entry increment the entrynumber and set partnumber to 0 00429 if ( partnumber == object->nParts - 1 || object->nParts == 0 ) 00430 { 00431 entrynumber++; 00432 entryindex++; 00433 partnumber = 0; 00434 SHPDestroyObject( object ); 00435 object=NULL; 00436 } 00437 else 00438 partnumber++; 00439 00440 if (nVertices==0) 00441 continue; 00442 #else 00443 // read in # points in segment 00444 if ( pdf_rdx( n_buff, (long) sizeof ( unsigned char ) * 2, in ) == 0 ) 00445 break; 00446 nVertices = ( n_buff[0] << 8 ) + n_buff[1]; 00447 if ( nVertices == 0 ) 00448 break; 00449 00450 pdf_rdx( buff, (long) sizeof ( unsigned char ) * 4 * nVertices, in ); 00451 if ( nVertices == 1 ) 00452 continue; 00453 00454 for ( j = i = 0; i < nVertices; i++, j += 2 ) 00455 { 00456 t = ( buff[j] << 8 ) + buff[j + 1]; 00457 bufx[i] = ( (PLFLT) t - OFFSET ) / SCALE; 00458 } 00459 for ( i = 0; i < nVertices; i++, j += 2 ) 00460 { 00461 t = ( buff[j] << 8 ) + buff[j + 1]; 00462 bufy[i] = ( (PLFLT) t - OFFSET ) / SCALE; 00463 } 00464 //set the min/max section lat/lon with extreme values 00465 //to be overwritten later 00466 minsectlon = 1000.; 00467 maxsectlon = -1000.; 00468 minsectlat = 1000.; 00469 maxsectlat = -1000.; 00470 00471 #endif 00472 00473 if( islatlon ) 00474 { 00475 //two obvious issues exist here with plotting longitudes: 00476 // 00477 //1) wraparound causing lines which go the wrong way round 00478 // the globe 00479 //2) some people plot lon from 0-360 deg, others from -180 - +180 00480 // 00481 //we can cure these problems by conditionally adding/subtracting 00482 //360 degrees to each data point in order to ensure that the 00483 //distance between adgacent points is always less than 180 00484 //degrees, then plotting up to 2 out of 5 copies of the data 00485 //each separated by 360 degrees. 00486 00487 //arrays of pointers to the starts of each section of data that 00488 //has been split due to longitude wrapping, and an array of ints 00489 //to hold their lengths. Start with splitx and splity having one 00490 //element pointing to the beginning of bufx and bufy 00491 splitx = ( PLFLT** ) malloc ( sizeof ( PLFLT* ) ); 00492 splity = ( PLFLT** ) malloc ( sizeof ( PLFLT* ) ); 00493 //lengths of the split sections 00494 splitsectionlengths = ( int* ) malloc ( sizeof ( size_t ) ); 00495 if(!splitx || ! splity || ! splitsectionlengths) 00496 { 00497 plabort("Could not allocate memory for longitudinally split map data"); 00498 free( filename ); 00499 free ( bufx ); 00500 free ( bufy ); 00501 free ( splitx ); 00502 free ( splity ); 00503 free ( splitsectionlengths ); 00504 return; 00505 } 00506 splitsectionlengths[0]=nVertices; 00507 nsplitsections=1; 00508 splitx[0]=bufx; 00509 splity[0]=bufy; 00510 00511 //set the min/max lats/lons 00512 minsectlon = MIN( minsectlon, bufx[0] ); 00513 maxsectlon = MAX( minsectlon, bufx[0] ); 00514 minsectlat = MIN( minsectlat, bufy[0] ); 00515 maxsectlat = MAX( maxsectlat, bufy[0] ); 00516 //ensure our lat and lon are on 0-360 grid and split the 00517 //data where it wraps. 00518 rebaselon( &bufx[0], ( minlong + maxlong ) / 2.0 ); 00519 for ( i = 1; i < nVertices; i++ ) 00520 { 00521 //put lon into 0-360 degree range 00522 rebaselon( &bufx[i], ( minlong + maxlong ) / 2.0 ); 00523 00524 //check if the previous point is more than 180 degrees away 00525 if ( bufx[i - 1] - bufx[i] > 180. || bufx[i - 1] - bufx[i] < -180.) 00526 { 00527 //check if the map transform deals with wrapping itself, e.g. in a polar projection 00528 //in this case give one point overlap to the sections so that lines are contiguous 00529 if( checkwrap( mapform, bufx[i], bufy[i] ) ) 00530 { 00531 appendresult += appendfltptr( &splitx, nsplitsections, bufx+i ); 00532 appendresult += appendfltptr( &splity, nsplitsections, bufy+i ); 00533 appendresult += appendint( &splitsectionlengths, nsplitsections, nVertices - i ); 00534 splitsectionlengths[nsplitsections-1] -= splitsectionlengths[nsplitsections] - 1; 00535 nsplitsections++; 00536 } 00537 //if the transform doesn't deal with wrapping then allow 2 points overlap to fill in the 00538 //edges 00539 else 00540 { 00541 appendresult += appendfltptr( &splitx, nsplitsections, bufx+i-1 ); 00542 appendresult += appendfltptr( &splity, nsplitsections, bufy+i-1 ); 00543 appendresult += appendint( &splitsectionlengths, nsplitsections, nVertices - i + 1 ); 00544 splitsectionlengths[nsplitsections-1] -= splitsectionlengths[nsplitsections] - 2; 00545 nsplitsections++; 00546 } 00547 if( appendresult > 0 ) 00548 { 00549 plabort("Could not allocate memory for appending to longitudinally split map data"); 00550 free( filename ); 00551 free ( bufx ); 00552 free ( bufy ); 00553 free ( splitx ); 00554 free ( splity ); 00555 free ( splitsectionlengths ); 00556 return; 00557 } 00558 00559 } 00560 00561 //update the mins and maxs 00562 minsectlon = MIN( minsectlon, bufx[i] ); 00563 maxsectlon = MAX( minsectlon, bufx[i] ); 00564 minsectlat = MIN( minsectlat, bufy[i] ); 00565 maxsectlat = MAX( maxsectlat, bufy[i] ); 00566 } 00567 00568 //check if the latitude and longitude range means we need to plot this section 00569 if ( ( maxsectlat > minlat ) && ( minsectlat < maxlat ) 00570 && ( maxsectlon > minlong ) && ( minsectlon < maxlong )) 00571 { 00572 //plot each split in turn, now is where we deal with the end points to 00573 //ensure we draw to the edge of the map 00574 for( i = 0; i < nsplitsections; ++i ) 00575 { 00576 //check if the first 2 or last 1 points of the split section need 00577 //wrapping and add or subtract 360 from them. Note that when the next 00578 //section is drawn the code below will undo this if needed 00579 if( splitsectionlengths[i] > 2 ) 00580 { 00581 if ( splitx[i][1] - splitx[i][2] > 180. ) 00582 splitx[i][1] -= 360.0; 00583 else if ( splitx[i][1] - splitx[i][2] < -180. ) 00584 splitx[i][1] += 360.0; 00585 } 00586 00587 if ( splitx[i][0] - splitx[i][1] > 180. ) 00588 splitx[i][0] -= 360.0; 00589 else if ( splitx[i][0] - splitx[i][1] < -180. ) 00590 splitx[i][0] += 360.0; 00591 00592 if ( splitx[i][splitsectionlengths[i] - 2] - splitx[i][splitsectionlengths[i] - 1] > 180. ) 00593 splitx[i][splitsectionlengths[i] - 1] += 360.0; 00594 else if ( splitx[i][splitsectionlengths[i] - 2] - splitx[i][splitsectionlengths[i] - 1] < -180. ) 00595 splitx[i][splitsectionlengths[i] - 1] -= 360.0; 00596 00597 //save the last 2 points - they will be needed by the next 00598 //split section and will be overwritten by the mapform 00599 lastsplitpointx = splitx[i][splitsectionlengths[i]-1]; 00600 lastsplitpointy = splity[i][splitsectionlengths[i]-1]; 00601 penultimatesplitpointx = splitx[i][splitsectionlengths[i]-2]; 00602 penultimatesplitpointy = splity[i][splitsectionlengths[i]-2]; 00603 00604 //draw the split section 00605 drawmapdata(mapform, shapetype, splitsectionlengths[i], splitx[i], splity[i], dx, dy, just, text ); 00606 00607 for(j=1; j<splitsectionlengths[i]; ++j) 00608 { 00609 if((splitx[i][j]<200.0 && splitx[i][j-1] > 260.0) || ( splitx[i][j-1]<200.0 && splitx[i][j] > 260.0 )) 00610 plwarn("wrapping error"); 00611 } 00612 00613 //restore the last 2 points 00614 splitx[i][splitsectionlengths[i]-1] = lastsplitpointx; 00615 splity[i][splitsectionlengths[i]-1] = lastsplitpointy; 00616 splitx[i][splitsectionlengths[i]-2] = penultimatesplitpointx; 00617 splity[i][splitsectionlengths[i]-2] = penultimatesplitpointy; 00618 } 00619 } 00620 } 00621 else 00622 { 00623 drawmapdata(mapform, shapetype, nVertices, bufx, bufy, dx, dy, just, text ); 00624 } 00625 00626 free( splitx ); 00627 free( splity ); 00628 free( splitsectionlengths ); 00629 00630 00631 } 00632 // Close map file 00633 #ifdef HAVE_SHAPELIB 00634 SHPClose( in ); 00635 #else 00636 pdf_close( in ); 00637 #endif 00638 00639 //free memory 00640 free( bufx ); 00641 free( bufy ); 00642 free( filename ); 00643 #else // defined (HAVE_SHAPELIB) || defined (PL_DEPRECATED) 00644 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" ); 00645 #endif // defined (HAVE_SHAPELIB) || defined (PL_DEPRECATED) 00646 } 00647 00648 00649 00650 //-------------------------------------------------------------------------- 00651 // void plmap(void (*mapform)(PLINT, PLFLT *, PLFLT *), const char *type, 00652 // PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat); 00653 // 00654 // plot continental outline in world coordinates 00655 // 00656 // v1.4: machine independant version 00657 // v1.3: replaced plcontinent by plmap, added plmeridians 00658 // v1.2: 2 arguments: mapform, type of plot 00659 // 00660 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the 00661 // coordinate longitudes and latitudes to a plot coordinate system. By 00662 // using this transform, we can change from a longitude, latitude 00663 // coordinate to a polar stereographic project, for example. Initially, 00664 // x[0]..[n-1] are the longitudes and y[0]..y[n-1] are the corresponding 00665 // latitudes. After the call to mapform(), x[] and y[] should be replaced 00666 // by the corresponding plot coordinates. If no transform is desired, 00667 // mapform can be replaced by NULL. 00668 // 00669 // type is a character string. The value of this parameter determines the 00670 // type of background. The possible values are, 00671 // 00672 // "globe" continental outlines 00673 // "usa" USA and state boundaries 00674 // "cglobe" continental outlines and countries 00675 // "usaglobe" USA, state boundaries and continental outlines 00676 // alternatively the filename of a shapefile can be passed if PLplot has 00677 // been compiled with shapelib. In this case either the base name of the 00678 // file can be passed or the filename including the .shp or .shx suffix. 00679 // Only the .shp and .shx files are used. 00680 // 00681 // minlong, maxlong are the values of the longitude on the left and right 00682 // side of the plot, respectively. The value of minlong must be less than 00683 // the values of maxlong, and the values of maxlong-minlong must be less 00684 // or equal to 360. 00685 // 00686 // minlat, maxlat are the minimum and maximum latitudes to be plotted on 00687 // the background. One can always use -90.0 and 90.0 as the boundary 00688 // outside the plot window will be automatically eliminated. However, the 00689 // program will be faster if one can reduce the size of the background 00690 // plotted. 00691 //-------------------------------------------------------------------------- 00692 00693 void 00694 plmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type, 00695 PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat ) 00696 { 00697 drawmap( mapform, type, 0.0, 0.0, SHPT_ARC, 0.0, NULL, minlong, maxlong, 00698 minlat, maxlat, NULL, 0 ); 00699 } 00700 00701 //-------------------------------------------------------------------------- 00702 // void plmapline( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 00703 // const char *type, PLFLT minlong, PLFLT maxlong, PLFLT minlat, 00704 // PLFLT maxlat, int* plotentries, int nplotentries); 00705 00706 //New version of plmap which allows us to specify which items in a shapefile 00707 //we want to use. parameters are as above but with the plotentries being an 00708 //array containing the indices of the elements we wish to draw and 00709 //nplotentries being the number of items in plotentries. 00710 //If shapefile access was not built into plplot then plotentries and 00711 //nplotentries are ignored. If plotentries is null than all entries are 00712 //drawn and nplotentries is ignored. 00713 //The name distiguishes it from other functions which plot points/text and 00714 //polygons, but note that the type of element in the shapefile need not 00715 //match the type of element drawn - i.e. arc elements from a shapefile could 00716 //be drawn as points using the plmaptex function. 00717 //-------------------------------------------------------------------------- 00718 void 00719 plmapline( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type, 00720 PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, 00721 int* plotentries, int nplotentries) 00722 { 00723 drawmap( mapform, type, 0.0, 0.0, SHPT_ARC, 0.0, "", minlong, maxlong, 00724 minlat, maxlat, plotentries, nplotentries ); 00725 } 00726 00727 //-------------------------------------------------------------------------- 00728 // void plmapstring( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 00729 // const char *type, PLFLT just, const char *string, 00730 // PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, 00731 // int* plotentries, int nplotentries); 00732 // 00733 //As per plmapline but plots symbols. The map equivalent of plstring. string 00734 //has the same meaning as in plstring. 00735 //------------------------------------------------------------------------- 00736 void 00737 plmapstring( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 00738 const char *type, const char *string, 00739 PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, 00740 int* plotentries, int nplotentries) 00741 { 00742 drawmap( mapform, type, 1.0, 0.0, SHPT_POINT, 0.5, string, minlong, maxlong, 00743 minlat, maxlat, plotentries, nplotentries ); 00744 } 00745 00746 //-------------------------------------------------------------------------- 00747 // void plmaptex( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 00748 // const char *type, PLFLT dx, PLFLT dy PLFLT just, const char *text, 00749 // PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, 00750 // int* plotentries, int nplotentries); 00751 // 00752 //As per plmapline but plots text. The map equivalent of plptex. dx, dy, 00753 //just and text have the same meaning as in plptex. 00754 //------------------------------------------------------------------------- 00755 void 00756 plmaptex( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 00757 const char *type, PLFLT dx, PLFLT dy, PLFLT just, const char *text, 00758 PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat, 00759 int plotentry) 00760 { 00761 drawmap( mapform, type, dx, dy, SHPT_POINT, just, text, minlong, maxlong, 00762 minlat, maxlat, &plotentry, 1 ); 00763 } 00764 00765 //-------------------------------------------------------------------------- 00766 // void plmapfill( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 00767 // const char *type, PLFLT minlong, PLFLT maxlong, PLFLT minlat, 00768 // PLFLT maxlat, int* plotentries, int nplotentries); 00769 // 00770 //As per plmapline but plots a filled polygon. The map equivalent to 00771 //plfill. Uses the pattern defined by plsty or plpat. 00772 //------------------------------------------------------------------------- 00773 void 00774 plmapfill( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 00775 const char *type, PLFLT minlong, PLFLT maxlong, PLFLT minlat, 00776 PLFLT maxlat, int* plotentries, int nplotentries) 00777 { 00778 drawmap( mapform, type, 0.0, 0.0, SHPT_POLYGON, 0.0, NULL, minlong, maxlong, 00779 minlat, maxlat, plotentries, nplotentries ); 00780 } 00781 00782 //-------------------------------------------------------------------------- 00783 // void plmeridians(void (*mapform)(PLINT, PLFLT *, PLFLT *), 00784 // PLFLT dlong, PLFLT dlat, PLFLT minlong, PLFLT maxlong, 00785 // PLFLT minlat, PLFLT maxlat); 00786 // 00787 // Plot the latitudes and longitudes on the background. The lines 00788 // are plotted in the current color and line style. 00789 // 00790 // mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the 00791 // coordinate longitudes and latitudes to a plot coordinate system. By 00792 // using this transform, we can change from a longitude, latitude 00793 // coordinate to a polar stereographic project, for example. Initially, 00794 // x[0]..x[n-1] are the longitudes and y[0]..y[n-1] are the corresponding 00795 // latitudes. After the call to mapform(), x[] and y[] should be replaced 00796 // by the corresponding plot coordinates. If no transform is desired, 00797 // mapform can be replaced by NULL. 00798 // 00799 // dlat, dlong are the interval in degrees that the latitude and longitude 00800 // lines are to be plotted. 00801 // 00802 // minlong, maxlong are the values of the longitude on the left and right 00803 // side of the plot, respectively. The value of minlong must be less than 00804 // the values of maxlong, and the values of maxlong-minlong must be less 00805 // or equal to 360. 00806 // 00807 // minlat, maxlat are the minimum and maximum latitudes to be plotted on 00808 // the background. One can always use -90.0 and 90.0 as the boundary 00809 // outside the plot window will be automatically eliminated. However, the 00810 // program will be faster if one can reduce the size of the background 00811 // plotted. 00812 //-------------------------------------------------------------------------- 00813 00814 #define NSEG 100 00815 00816 void 00817 plmeridians( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), 00818 PLFLT dlong, PLFLT dlat, 00819 PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat ) 00820 { 00821 PLFLT yy, xx, temp, x[2], y[2], dx, dy; 00822 00823 if ( minlong > maxlong ) 00824 { 00825 temp = minlong; 00826 minlong = maxlong; 00827 maxlong = temp; 00828 } 00829 if ( minlat > maxlat ) 00830 { 00831 temp = minlat; 00832 minlat = maxlat; 00833 maxlat = temp; 00834 } 00835 dx = ( maxlong - minlong ) / NSEG; 00836 dy = ( maxlat - minlat ) / NSEG; 00837 00838 // latitudes 00839 00840 for ( yy = dlat * ceil( minlat / dlat ); yy <= maxlat; yy += dlat ) 00841 { 00842 if ( mapform == NULL ) 00843 { 00844 plpath( NSEG, minlong, yy, maxlong, yy ); 00845 } 00846 else 00847 { 00848 for ( xx = minlong; xx < maxlong; xx += dx ) 00849 { 00850 y[0] = y[1] = yy; 00851 x[0] = xx; 00852 x[1] = xx + dx; 00853 ( *mapform )( 2, x, y ); 00854 plline( 2, x, y ); 00855 } 00856 } 00857 } 00858 00859 // longitudes 00860 00861 for ( xx = dlong * ceil( minlong / dlong ); xx <= maxlong; xx += dlong ) 00862 { 00863 if ( mapform == NULL ) 00864 { 00865 plpath( NSEG, xx, minlat, xx, maxlat ); 00866 } 00867 else 00868 { 00869 for ( yy = minlat; yy < maxlat; yy += dy ) 00870 { 00871 x[0] = x[1] = xx; 00872 y[0] = yy; 00873 y[1] = yy + dy; 00874 ( *mapform )( 2, x, y ); 00875 plline( 2, x, y ); 00876 } 00877 } 00878 } 00879 } 00880 00881 //-------------------------------------------------------------------------- 00882 // SHPHandle OpenShapeFile(fn) 00883 // 00897 //-------------------------------------------------------------------------- 00898 #ifdef HAVE_SHAPELIB 00899 #ifdef HAVE_SAHOOKS 00900 // Our thanks to Frank Warmerdam, the developer of shapelib for suggesting 00901 // this approach for quieting shapelib "Unable to open" error messages. 00902 static 00903 void CustomErrors( const char *message ) 00904 { 00905 if ( strstr( message, "Unable to open" ) == NULL ) 00906 fprintf( stderr, "%s\n", message ); 00907 } 00908 #endif 00909 00910 SHPHandle 00911 OpenShapeFile( const char *fn ) 00912 { 00913 SHPHandle file; 00914 char *fs = NULL, *dn = NULL; 00915 #ifdef HAVE_SAHOOKS 00916 SAHooks sHooks; 00917 00918 SASetupDefaultHooks( &sHooks ); 00919 sHooks.Error = CustomErrors; 00920 #else 00921 // Using ancient version of shapelib without SAHooks or SHPOpenLL. 00922 // For this case live with the misleading "Unable to open" error 00923 // messages. 00924 // int sHooks; 00925 #define SHPOpenLL( a, b, c ) SHPOpen( a, b ) 00926 #endif 00927 00928 //*** search build tree *** 00929 00930 if ( plInBuildTree() == 1 ) 00931 { 00932 plGetName( SOURCE_DIR, "data", fn, &fs ); 00933 00934 if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL ) 00935 goto done; 00936 } 00937 00938 //*** search PLPLOT_LIB_ENV = $(PLPLOT_LIB) *** 00939 00940 #if defined ( PLPLOT_LIB_ENV ) 00941 if ( ( dn = getenv( PLPLOT_LIB_ENV ) ) != NULL ) 00942 { 00943 plGetName( dn, "", fn, &fs ); 00944 00945 if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL ) 00946 goto done; 00947 fprintf( stderr, PLPLOT_LIB_ENV "=\"%s\"\n", dn ); // what IS set? 00948 } 00949 #endif // PLPLOT_LIB_ENV 00950 00951 //*** search current directory *** 00952 00953 if ( ( file = SHPOpenLL( fn, "rb", &sHooks ) ) != NULL ) 00954 { 00955 pldebug( "OpenShapeFile", "Found file %s in current directory.\n", fn ); 00956 free_mem( fs ); 00957 return ( file ); 00958 } 00959 00960 //*** search PLPLOT_HOME_ENV/lib = $(PLPLOT_HOME)/lib *** 00961 00962 #if defined ( PLPLOT_HOME_ENV ) 00963 if ( ( dn = getenv( PLPLOT_HOME_ENV ) ) != NULL ) 00964 { 00965 plGetName( dn, "lib", fn, &fs ); 00966 00967 if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL ) 00968 goto done; 00969 fprintf( stderr, PLPLOT_HOME_ENV "=\"%s\"\n", dn ); // what IS set? 00970 } 00971 #endif // PLPLOT_HOME_ENV/lib 00972 00973 //*** search installed location *** 00974 00975 #if defined ( DATA_DIR ) 00976 plGetName( DATA_DIR, "", fn, &fs ); 00977 00978 if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL ) 00979 goto done; 00980 #endif // DATA_DIR 00981 00982 //*** search hardwired location *** 00983 00984 #ifdef PLLIBDEV 00985 plGetName( PLLIBDEV, "", fn, &fs ); 00986 00987 if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL ) 00988 goto done; 00989 #endif // PLLIBDEV 00990 00991 //*** not found, give up *** 00992 pldebug( "OpenShapeFile", "File %s not found.\n", fn ); 00993 free_mem( fs ); 00994 return NULL; 00995 00996 done: 00997 pldebug( "OpenShapeFile", "SHPOpen successfully opened two files with basename %s\n", fs ); 00998 free_mem( fs ); 00999 return ( file ); 01000 } 01001 #endif