PLplot 5.9.6
plmap.c
00001 /* $Id$
00002  *
00003  *      Continental Outline and Political Boundary Backgrounds
00004  *
00005  *      Some plots need a geographical background such as the global
00006  *      surface temperatures or the population density.  The routine
00007  *      plmap() will draw one of the following backgrounds: continental
00008  *      outlines, political boundaries, the United States, and the United
00009  *      States with the continental outlines.  The routine plmeridians()
00010  *      will add the latitudes and longitudes to the background.  After
00011  *      the background has been drawn, one can use a contour routine or a
00012  *      symbol plotter to finish off the plot.
00013  *
00014  *      Copyright (C) 1991, 1993, 1994  Wesley Ebisuzaki
00015  *      Copyright (C) 1994, 2000, 2001  Maurice LeBrun
00016  *      Copyright (C) 1999  Geoffrey Furnish
00017  *      Copyright (C) 2000, 2001, 2002  Alan W. Irwin
00018  *      Copyright (C) 2001  Andrew Roach
00019  *      Copyright (C) 2001, 2004  Rafael Laboissiere
00020  *      Copyright (C) 2002  Vincent Darley
00021  *      Copyright (C) 2003  Joao Cardoso
00022  *
00023  *      This file is part of PLplot.
00024  *
00025  *      PLplot is free software; you can redistribute it and/or modify
00026  *      it under the terms of the GNU Library General Public License
00027  *      as published by the Free Software Foundation; version 2 of the
00028  *      License.
00029  *
00030  *      PLplot is distributed in the hope that it will be useful, but
00031  *      WITHOUT ANY WARRANTY; without even the implied warranty of
00032  *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00033  *      GNU Library General Public License for more details.
00034  *
00035  *      You should have received a copy of the GNU Library General Public
00036  *      License along with this library; if not, write to the Free Software
00037  *      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301
00038  *      USA
00039  *
00040  */
00041 
00042 #include "plplotP.h"
00043 
00044 /*----------------------------------------------------------------------*\
00045  * void plmap(void (*mapform)(PLINT, PLFLT *, PLFLT *), const char *type,
00046  *            PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat);
00047  *
00048  * plot continental outline in world coordinates
00049  *
00050  * v1.4: machine independant version
00051  * v1.3: replaced plcontinent by plmap, added plmeridians
00052  * v1.2: 2 arguments:  mapform, type of plot
00053  *
00054  * mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
00055  * coordinate longitudes and latitudes to a plot coordinate system.  By
00056  * using this transform, we can change from a longitude, latitude
00057  * coordinate to a polar stereographic project, for example.  Initially,
00058  * x[0]..[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
00059  * latitudes.  After the call to mapform(), x[] and y[] should be replaced
00060  * by the corresponding plot coordinates.  If no transform is desired,
00061  * mapform can be replaced by NULL.
00062  *
00063  * type is a character string. The value of this parameter determines the
00064  * type of background. The possible values are,
00065  *
00066  *      "globe"     continental outlines
00067  *      "usa"       USA and state boundaries
00068  *      "cglobe"    continental outlines and countries
00069  *      "usaglobe"  USA, state boundaries and continental outlines
00070  *
00071  * minlong, maxlong are the values of the longitude on the left and right
00072  * side of the plot, respectively. The value of minlong must be less than
00073  * the values of maxlong, and the values of maxlong-minlong must be less
00074  * or equal to 360.
00075  *
00076  * minlat, maxlat are the minimum and maximum latitudes to be plotted on
00077  * the background.  One can always use -90.0 and 90.0 as the boundary
00078  * outside the plot window will be automatically eliminated.  However, the
00079  * program will be faster if one can reduce the size of the background
00080  * plotted.
00081  \*----------------------------------------------------------------------*/
00082 
00083 #define MAP_FILE    ".map"
00084 #define OFFSET      ( 180 * 100 )
00085 #define SCALE       100.0
00086 #define W_BUFSIZ    ( 32 * 1024 )
00087 
00088 void
00089 plmap( void ( *mapform )( PLINT, PLFLT *, PLFLT * ), const char *type,
00090        PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
00091 {
00092     PLINT            wrap, sign;
00093     int              i, j;
00094     PLFLT            bufx[200], bufy[200], x[2], y[2];
00095     short int        test[400];
00096     register PDFstrm *in;
00097     char             filename[100];
00098 
00099     unsigned char    n_buff[2], buff[800];
00100     int              n;
00101     long int         t;
00102 
00103     /*
00104      * read map outline
00105      */
00106     strncpy( filename, type, 99 );
00107     filename[99] = '\0';
00108     strncat( filename, MAP_FILE, 99 - strlen( filename ) );
00109 
00110     if ( ( in = plLibOpenPdfstrm( filename ) ) == NULL )
00111         return;
00112 
00113     for (;; )
00114     {
00115         /* read in # points in segment */
00116         if ( pdf_rdx( n_buff, sizeof ( unsigned char ) * 2, in ) == 0 )
00117             break;
00118         n = ( n_buff[0] << 8 ) + n_buff[1];
00119         if ( n == 0 )
00120             break;
00121 
00122         pdf_rdx( buff, sizeof ( unsigned char ) * 4 * n, in );
00123         if ( n == 1 )
00124             continue;
00125 
00126         for ( j = i = 0; i < n; i++, j += 2 )
00127         {
00128             t       = ( buff[j] << 8 ) + buff[j + 1];
00129             bufx[i] = ( t - OFFSET ) / SCALE;
00130         }
00131         for ( i = 0; i < n; i++, j += 2 )
00132         {
00133             t       = ( buff[j] << 8 ) + buff[j + 1];
00134             bufy[i] = ( t - OFFSET ) / SCALE;
00135         }
00136 
00137         for ( i = 0; i < n; i++ )
00138         {
00139             while ( bufx[i] < minlong )
00140             {
00141                 bufx[i] += 360.0;
00142             }
00143             while ( bufx[i] > maxlong )
00144             {
00145                 bufx[i] -= 360.0;
00146             }
00147         }
00148 
00149         /* remove last 2 points if both outside of domain and won't plot */
00150 
00151 /* AR: 18/11/01
00152  *       I have commented out the next section which supposedly
00153  *       removes points that do not plot within the domain.
00154  *
00155  *       This code appears at any rate to be superseded by the next
00156  *       block of code that checks for wrapping problems. Removing
00157  *       this code seems to have fixed up the problems with mapping
00158  *       function, but I do not wish to delete it outright just now in
00159  *       case I have managed to miss something.
00160  */
00161 
00162 /*  while (n > 1) {
00163  *          if ((bufx[n-1] < minlong && bufx[n-2] < minlong) ||
00164  *              (bufx[n-1] > maxlong && bufx[n-2] > maxlong) ||
00165  *              (bufy[n-1] < minlat && bufy[n-2] < minlat) ||
00166  *              (bufy[n-1] > maxlat && bufy[n-2] > maxlat)) {
00167  *              n--;
00168  *          }
00169  *          else {
00170  *              break;
00171  *          }
00172  *      }
00173  *      if (n <= 1) continue;
00174  */
00175 
00176         if ( mapform != NULL )
00177             ( *mapform )( n, bufx, bufy );                    /* moved transformation to here   */
00178                                                               /* so bound-checking worked right */
00179 
00180         wrap = 0;
00181         /* check for wrap around problems */
00182         for ( i = 0; i < n - 1; i++ )
00183         {
00184             /* jc: this code is wrong, as the bufx/y are floats that are
00185              * converted to ints before abs() is called. Thus, small
00186              * differences are masked off. The proof that it is wrong is
00187              * that when replacing abs() with fabs(), as it should be,
00188              * the code works the wrong way. What a proof :-)
00189              *
00190              * test[i] = abs((bufx[i]-bufx[i+1])) > abs(bufy[i]/3);
00191              *
00192              * If the intended behaviour is *really* that, than an
00193              * explicit cast should be used to help other programmers, as
00194              * is done bellow!!!
00195              */
00196 
00197             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" */
00198             if ( test[i] )
00199                 wrap = 1;
00200         }
00201 
00202         if ( wrap == 0 )
00203         {
00204             plline( n, bufx, bufy );
00205         }
00206         else
00207         {
00208             for ( i = 0; i < n - 1; i++ )
00209             {
00210                 x[0] = bufx[i];
00211                 x[1] = bufx[i + 1];
00212                 y[0] = bufy[i];
00213                 y[1] = bufy[i + 1];
00214                 if ( test[i] == 0 )
00215                 {
00216                     plline( 2, x, y );
00217                 }
00218                 else    /* this code seems to supercede the block commented out above */
00219                 {
00220                     /* segment goes off the edge */
00221                     sign  = ( x[1] > x[0] ) ? 1 : -1;
00222                     x[1] -= sign * 360.0;
00223                     plline( 2, x, y );
00224                     x[0]  = bufx[i];
00225                     x[1]  = bufx[i + 1];
00226                     y[0]  = bufy[i];
00227                     y[1]  = bufy[i + 1];
00228                     x[0] += sign * 360.0;
00229                     plline( 2, x, y );
00230                 }
00231             }
00232         }
00233     }
00234     /* Close map file */
00235     pdf_close( in );
00236 }
00237 
00238 /*----------------------------------------------------------------------*\
00239  * void plmeridians(void (*mapform)(PLINT, PLFLT *, PLFLT *),
00240  *          PLFLT dlong, PLFLT dlat, PLFLT minlong, PLFLT maxlong,
00241  *          PLFLT minlat, PLFLT maxlat);
00242  *
00243  * Plot the latitudes and longitudes on the background.  The lines
00244  * are plotted in the current color and line style.
00245  *
00246  * mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
00247  * coordinate longitudes and latitudes to a plot coordinate system.  By
00248  * using this transform, we can change from a longitude, latitude
00249  * coordinate to a polar stereographic project, for example.  Initially,
00250  * x[0]..x[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
00251  * latitudes.  After the call to mapform(), x[] and y[] should be replaced
00252  * by the corresponding plot coordinates.  If no transform is desired,
00253  * mapform can be replaced by NULL.
00254  *
00255  * dlat, dlong are the interval in degrees that the latitude and longitude
00256  * lines are to be plotted.
00257  *
00258  * minlong, maxlong are the values of the longitude on the left and right
00259  * side of the plot, respectively. The value of minlong must be less than
00260  * the values of maxlong, and the values of maxlong-minlong must be less
00261  * or equal to 360.
00262  *
00263  * minlat, maxlat are the minimum and maximum latitudes to be plotted on
00264  * the background.  One can always use -90.0 and 90.0 as the boundary
00265  * outside the plot window will be automatically eliminated.  However, the
00266  * program will be faster if one can reduce the size of the background
00267  * plotted.
00268  \*----------------------------------------------------------------------*/
00269 
00270 #define NSEG    100
00271 
00272 void
00273 plmeridians( void ( *mapform )( PLINT, PLFLT *, PLFLT * ),
00274              PLFLT dlong, PLFLT dlat,
00275              PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
00276 {
00277     PLFLT yy, xx, temp, x[2], y[2], dx, dy;
00278 
00279     if ( minlong > maxlong )
00280     {
00281         temp    = minlong;
00282         minlong = maxlong;
00283         maxlong = temp;
00284     }
00285     if ( minlat > maxlat )
00286     {
00287         temp   = minlat;
00288         minlat = maxlat;
00289         maxlat = temp;
00290     }
00291     dx = ( maxlong - minlong ) / NSEG;
00292     dy = ( maxlat - minlat ) / NSEG;
00293 
00294     /* latitudes */
00295 
00296     for ( yy = dlat * ceil( minlat / dlat ); yy <= maxlat; yy += dlat )
00297     {
00298         if ( mapform == NULL )
00299         {
00300             plpath( NSEG, minlong, yy, maxlong, yy );
00301         }
00302         else
00303         {
00304             for ( xx = minlong; xx < maxlong; xx += dx )
00305             {
00306                 y[0] = y[1] = yy;
00307                 x[0] = xx;
00308                 x[1] = xx + dx;
00309                 ( *mapform )( 2, x, y );
00310                 plline( 2, x, y );
00311             }
00312         }
00313     }
00314 
00315     /* longitudes */
00316 
00317     for ( xx = dlong * ceil( minlong / dlong ); xx <= maxlong; xx += dlong )
00318     {
00319         if ( mapform == NULL )
00320         {
00321             plpath( NSEG, xx, minlat, xx, maxlat );
00322         }
00323         else
00324         {
00325             for ( yy = minlat; yy < maxlat; yy += dy )
00326             {
00327                 x[0] = x[1] = xx;
00328                 y[0] = yy;
00329                 y[1] = yy + dy;
00330                 ( *mapform )( 2, x, y );
00331                 plline( 2, x, y );
00332             }
00333         }
00334     }
00335 }
 All Data Structures Files Functions