Index: src/plmap.c =================================================================== --- src/plmap.c (revision 12241) +++ src/plmap.c (working copy) @@ -39,8 +39,22 @@ // // +#define NEED_PLDEBUG 1 + #include "plplotP.h" + +// Used by any external init code to suggest a path +static char PLDLLIMPEXP * plplotLibDir; + +#ifdef HAVE_SHAPELIB +#include + +SHPHandle +OpenShapeFile( const char *fn); + +#endif + //-------------------------------------------------------------------------- // void plmap(void (*mapform)(PLINT, PLFLT *, PLFLT *), const char *type, // PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat); @@ -80,25 +94,52 @@ // plotted. //-------------------------------------------------------------------------- +#ifdef HAVE_SHAPELIB +#define MAP_FILE "" +#define OpenMap OpenShapeFile +#define CloseMap SHPClose +#else #define MAP_FILE ".map" +#define OpenMap plLibOpenPdfstrm +#define CloseMap pdf_close +#endif + #define OFFSET ( 180 * 100 ) #define SCALE 100.0 #define W_BUFSIZ ( 32 * 1024 ) +#define NCOPIES 5 //must be odd - original plus copies shifted by multiples of +/- 360 void plmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type, - PLFLT minlong, PLFLT maxlong, PLFLT PL_UNUSED( minlat ), PLFLT PL_UNUSED( maxlat ) ) + PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat ) { - PLINT wrap, sign; int i, j; - PLFLT bufx[200], bufy[200], x[2], y[2]; - short int test[400]; - register PDFstrm *in; char filename[100]; + int n; + PLFLT minsectlon, maxsectlon, minsectlat, maxsectlat; + const int mid=NCOPIES/2+1; +#ifdef HAVE_SHAPELIB + SHPHandle in; + int nentries; + int nparts; + int entrynumber=0; + int partnumber=0; + int shapetype; + double mins[4]; + double maxs[4]; + SHPObject *object=NULL; + double *bufxraw; + double *bufyraw; + PLFLT **bufx=NULL, **bufy=NULL; + int bufsize=0; +#else + register PDFstrm *in; + PLFLT bufx[NCOPIES][200], bufy[NCOPIES][200]; unsigned char n_buff[2], buff[800]; - int n; long int t; + int k; +#endif // // read map outline @@ -107,11 +148,79 @@ filename[99] = '\0'; strncat( filename, MAP_FILE, 99 - strlen( filename ) ); - if ( ( in = plLibOpenPdfstrm( filename ) ) == NULL ) - return; + if( ( in = OpenMap( filename ) ) == NULL ) + return; +#ifdef HAVE_SHAPELIB + SHPGetInfo( in, &nentries, &shapetype, mins, maxs ); + bufx = malloc( NCOPIES * sizeof(PLFLT*) ); + bufy = malloc( NCOPIES * sizeof(PLFLT*) ); + for( i = 0; i < NCOPIES; i++ ) + { + bufx[i] = NULL; + bufy[i] = NULL; + } +#endif for (;; ) { +#ifdef HAVE_SHAPELIB + //break condition if we've reached the end of the file + if ( entrynumber == nentries ) + break; + //if partnumber == 0 then we need to load the next object + if ( partnumber == 0 ) + { + object = SHPReadObject( in, entrynumber ); + nparts = object->nParts; + } + + //work out how many points are in the current part + if ( partnumber == ( nparts - 1 ) ) + n = object->nVertices - object->panPartStart[partnumber]; + else + n = object->panPartStart[partnumber+1] - object->panPartStart[partnumber]; + //allocate memory for the data + if ( n > bufsize ) + { + bufsize = n; + for ( i = 0; i < NCOPIES; i++) + { + if ( bufx[i] ) + free( bufx[i] ); + if ( bufy[i] ) + free( bufy[i] ); + bufx[i] = malloc( bufsize * sizeof(double) ); + bufy[i] = malloc( bufsize * sizeof(double) ); + } + } + + //point the plot buffer to the correct starting vertex + //and copy it to the PLFLT arrays + bufxraw = object->padfX + object->panPartStart[partnumber]; + bufyraw = object->padfY + object->panPartStart[partnumber]; + for ( i = 0; i < n; i++ ) + { + bufx[mid][i] = (PLFLT) bufxraw[i]; + for ( j = 0; j < NCOPIES; j++ ) + bufy[j][i] = (PLFLT) bufyraw[i]; + } + + //set the minlat/lon of the object + minsectlon = object->dfXMin; + maxsectlon = object->dfXMax; + minsectlat = object->dfYMin; + maxsectlat = object->dfYMax; + + //increment the partnumber or if we've reached the end of + //an entry increment the entrynumber and set partnumber to 0 + if ( partnumber == nparts -1 ) + { + entrynumber++; + partnumber = 0; + } + else + partnumber++; +#else // read in # points in segment if ( pdf_rdx( n_buff, (long) sizeof ( unsigned char ) * 2, in ) == 0 ) break; @@ -126,113 +235,97 @@ for ( j = i = 0; i < n; i++, j += 2 ) { t = ( buff[j] << 8 ) + buff[j + 1]; - bufx[i] = ( (PLFLT) t - OFFSET ) / SCALE; + bufx[mid][i] = ( (PLFLT) t - OFFSET ) / SCALE; } for ( i = 0; i < n; i++, j += 2 ) { t = ( buff[j] << 8 ) + buff[j + 1]; - bufy[i] = ( (PLFLT) t - OFFSET ) / SCALE; + bufy[0][i] = ( (PLFLT) t - OFFSET ) / SCALE; + for ( k = 1; k < NCOPIES; k++ ) + bufy[k][i] = bufy[0][i]; } - for ( i = 0; i < n; i++ ) - { - while ( bufx[i] < minlong ) - { - bufx[i] += 360.0; - } - while ( bufx[i] > maxlong ) - { - bufx[i] -= 360.0; - } - } + //set the min/max section lat/lon with extreme values + //to be overwritten later + minsectlon = 1000.; + maxsectlon = -1000.; + minsectlat = 1000.; + maxsectlat = -1000.; - // remove last 2 points if both outside of domain and won't plot - -// AR: 18/11/01 -// I have commented out the next section which supposedly -// removes points that do not plot within the domain. -// -// This code appears at any rate to be superseded by the next -// block of code that checks for wrapping problems. Removing -// this code seems to have fixed up the problems with mapping -// function, but I do not wish to delete it outright just now in -// case I have managed to miss something. -// - -// while (n > 1) { -// if ((bufx[n-1] < minlong && bufx[n-2] < minlong) || -// (bufx[n-1] > maxlong && bufx[n-2] > maxlong) || -// (bufy[n-1] < minlat && bufy[n-2] < minlat) || -// (bufy[n-1] > maxlat && bufy[n-2] > maxlat)) { -// n--; -// } -// else { -// break; -// } -// } -// if (n <= 1) continue; -// - - if ( mapform != NULL ) - ( *mapform )( n, bufx, bufy ); // moved transformation to here - // so bound-checking worked right - - wrap = 0; - // check for wrap around problems +#endif + //two obvious issues exist here with plotting longitudes: + // + //1) wraparound causing lines which go the wrong way round + // the globe + //2) some people plot lon from 0-360 deg, others from -180 - +180 + // + //we can cure these problems by conditionally adding/subtracting + //360 degrees to each data point in order to ensure that the + //distance between adgacent points is always less than 180 + //degrees, then plotting up to 2 out of 5 copies of the data + //each separated by 360 degrees. + for ( i = 0; i < n - 1; i++ ) { - // jc: this code is wrong, as the bufx/y are floats that are - // converted to ints before abs() is called. Thus, small - // differences are masked off. The proof that it is wrong is - // that when replacing abs() with fabs(), as it should be, - // the code works the wrong way. What a proof :-) - // - // test[i] = abs((bufx[i]-bufx[i+1])) > abs(bufy[i]/3); - // - // If the intended behaviour is *really* that, than an - // explicit cast should be used to help other programmers, as - // is done bellow!!! - // + if ( bufx[mid][i]-bufx[mid][i+1] > 180. ) + bufx[mid][i+1]+=360.; + else if ( bufx[mid][i]-bufx[mid][i+1] < -180. ) + bufx[mid][i+1]-=360.; + } + for ( i = 0; i < n; i++ ) + { + for ( j = 0; j < mid; j++ ) + bufx[j][i]=bufx[mid][i] + 360. * (PLFLT)( j - mid ); + for ( j = mid + 1; j < NCOPIES; j++ ) + bufx[j][i]=bufx[mid][i] + 360. * (PLFLT)( j - mid ); +#ifndef HAVE_SHAPELIB + minsectlon=MIN( minsectlon, bufx[mid][i] ); + maxsectlon=MAX( minsectlon, bufx[mid][i] ); + minsectlat=MIN( minsectlat, bufy[mid][i] ); + maxsectlat=MAX( minsectlat, bufy[mid][i] ); +#endif + } - test[i] = abs( (int) ( bufx[i] - bufx[i + 1] ) ) > abs( (int) bufy[i] / 3 ); // Changed this from 30 degrees so it is now "polar sensitive" - if ( test[i] ) - wrap = 1; - } + //check if the latitude range means we need to plot this section + if ( ( maxsectlat > minlat ) && ( minsectlat < maxlat ) ) + { + //check which of the translated maps fall within the + //range and transform and plot them - note more than one + //map may be needed due to wrapping + for (j = 0; j < NCOPIES; j++ ) + { + if( ( minsectlon + 360. * (PLFLT)( j - mid ) < maxlong ) + && ( maxsectlon + 360. * (PLFLT)( j - mid ) > minlong ) ) + { + if ( mapform != NULL ) + ( *mapform )( n, bufx[j], bufy[j] ); + plline( n, bufx[j], bufy[j] ); + } + } + } - if ( wrap == 0 ) - { - plline( n, bufx, bufy ); - } - else - { - for ( i = 0; i < n - 1; i++ ) - { - x[0] = bufx[i]; - x[1] = bufx[i + 1]; - y[0] = bufy[i]; - y[1] = bufy[i + 1]; - if ( test[i] == 0 ) - { - plline( 2, x, y ); - } - else // this code seems to supercede the block commented out above - { - // segment goes off the edge - sign = ( x[1] > x[0] ) ? 1 : -1; - x[1] -= sign * 360.0; - plline( 2, x, y ); - x[0] = bufx[i]; - x[1] = bufx[i + 1]; - y[0] = bufy[i]; - y[1] = bufy[i + 1]; - x[0] += sign * 360.0; - plline( 2, x, y ); - } - } - } + + +#ifdef HAVE_SHAPELIB + if ( partnumber == 0 ) + SHPDestroyObject(object); +#endif } // Close map file - pdf_close( in ); + CloseMap( in ); + +#ifdef HAVE_SHAPELIB + //free memory + for ( i = 0; i < NCOPIES; i++ ) + { + if ( bufx[i] ) + free( bufx[i] ); + if ( bufy[i] ) + free( bufy[i] ); + } + free( bufx ); + free( bufy ); +#endif } //-------------------------------------------------------------------------- @@ -333,3 +426,115 @@ } } } + +//-------------------------------------------------------------------------- +// SHPHandle OpenShapeFile(fn) +// +//! Returns a handle to a shapefile given filename fn. Content based on +//! plLibOpenPdfstrm in plctrl.c +//! Locations checked: +//! PLPLOT_LIB_ENV = $(PLPLOT_LIB) +//! current directory +//! PLPLOT_HOME_ENV/lib = $(PLPLOT_HOME)/lib +//! DATA_DIR +//! PLLIBDEV +//! +//! @param fn Name of the file. +//! +//! @Return handle to a shapefile on success or NULL if the file cannot be +//! found +//-------------------------------------------------------------------------- +#ifdef HAVE_SHAPELIB +SHPHandle +OpenShapeFile( const char *fn) +{ + SHPHandle file; + char *fs = NULL, *dn = NULL; + +//*** search build tree *** + + if ( plInBuildTree() == 1 ) + { + plGetName( SOURCE_DIR, "data", fn, &fs ); + + if ( ( file = SHPOpen( fs, "rb" ) ) != NULL ) + goto done; + } + +//*** search PLPLOT_LIB_ENV = $(PLPLOT_LIB) *** + +#if defined ( PLPLOT_LIB_ENV ) + if ( ( dn = getenv( PLPLOT_LIB_ENV ) ) != NULL ) + { + plGetName( dn, "", fn, &fs ); + + if ( ( file = SHPOpen( fs, "rb" ) ) != NULL ) + goto done; + fprintf( stderr, PLPLOT_LIB_ENV "=\"%s\"\n", dn ); // what IS set? + } +#endif // PLPLOT_LIB_ENV + +//*** search current directory *** + + if ( ( file = SHPOpen( fn, "rb" ) ) != NULL ) + { + pldebug( "plLibOpenPdfstr", "Found file %s in current directory.\n", fn ); + free_mem( fs ); + return ( file ); + } + +//*** search PLPLOT_HOME_ENV/lib = $(PLPLOT_HOME)/lib *** + +#if defined ( PLPLOT_HOME_ENV ) + if ( ( dn = getenv( PLPLOT_HOME_ENV ) ) != NULL ) + { + plGetName( dn, "lib", fn, &fs ); + + if ( ( file = SHPOpen( fs, "rb" ) ) != NULL ) + goto done; + fprintf( stderr, PLPLOT_HOME_ENV "=\"%s\"\n", dn ); // what IS set? + } +#endif // PLPLOT_HOME_ENV/lib + +//*** search installed location *** + +#if defined ( DATA_DIR ) + plGetName( DATA_DIR, "", fn, &fs ); + + if ( ( file = SHPOpen( fs, "rb" ) ) != NULL ) + goto done; +#endif // DATA_DIR + +//*** search hardwired location *** + +#ifdef PLLIBDEV + plGetName( PLLIBDEV, "", fn, &fs ); + + if ( ( file = pdf_fopen( fs, "rb" ) ) != NULL ) + goto done; +#endif // PLLIBDEV + +#ifdef macintosh + file = plMacLibOpen( fn ); + if ( file != NULL ) + goto done; +#endif // macintosh + + if ( plplotLibDir != NULL ) + { + plGetName( plplotLibDir, "", fn, &fs ); + if ( ( file = SHPOpen( fs, "rb" ) ) != NULL ) + goto done; + } + +//*** not found, give up *** + pldebug( "plLibOpenPdfstr", "File %s not found.\n", fn ); + free_mem( fs ); + return NULL; + +done: + pldebug( "plLibOpenPdfstr", "Found file %s\n", fs ); + free_mem( fs ); + return ( file ); +} +#endif \ No newline at end of file