Back to index

plt-scheme  4.2.1
plmap.c
Go to the documentation of this file.
00001 /* $Id: plmap.c,v 1.1 2004/03/01 20:54:52 cozmic Exp $
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  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., 59 Temple Place, Suite 330, Boston, MA  02111-1307
00038         USA
00039 
00040 */
00041 
00042 #include "plplotP.h"
00043 
00044 /*----------------------------------------------------------------------*\
00045  * void plmap(void (*mapform)(PLINT, PLFLT *, PLFLT *), 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 *), 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     strcpy(filename,type);
00107     strcat(filename,MAP_FILE);
00108 
00109     if ((in = plLibOpenPdfstrm(filename)) == NULL)
00110        return;
00111 
00112     for (;;) {
00113        /* read in # points in segment */
00114        if (pdf_rdx(n_buff, sizeof (unsigned char)* 2, in) == 0) break;
00115        n = (n_buff[0] << 8) + n_buff[1];
00116        if (n == 0) break;
00117 
00118        pdf_rdx(buff, sizeof (unsigned char)*4*n, in);
00119        if (n == 1) continue;
00120 
00121        for (j = i = 0; i < n; i++, j += 2) {
00122            t = (buff[j] << 8) + buff[j+1];
00123            bufx[i] = (t - OFFSET) / SCALE;
00124        }
00125        for (i = 0; i < n; i++, j += 2) {
00126            t = (buff[j] << 8) + buff[j+1];
00127            bufy[i] = (t - OFFSET) / SCALE;
00128        }
00129 
00130        for (i = 0; i < n; i++) {
00131            while (bufx[i] < minlong) {
00132               bufx[i] += 360.0;
00133            }
00134            while (bufx[i] > maxlong) {
00135               bufx[i] -= 360.0;
00136            }
00137        }
00138 
00139        /* remove last 2 points if both outside of domain and won't plot */
00140 
00141 /* AR: 18/11/01 
00142 *       I have commented out the next section which supposedly
00143 *       removes points that do not plot within the domain. 
00144 *       
00145 *       This code appears at any rate to be superseded by the next
00146 *       block of code that checks for wrapping problems. Removing
00147 *       this code seems to have fixed up the problems with mapping
00148 *       function, but I do not wish to delete it outright just now in
00149 *       case I have managed to miss something.
00150 */
00151 
00152 /*     while (n > 1) {
00153            if ((bufx[n-1] < minlong && bufx[n-2] < minlong) ||
00154               (bufx[n-1] > maxlong && bufx[n-2] > maxlong) ||
00155               (bufy[n-1] < minlat && bufy[n-2] < minlat) ||
00156               (bufy[n-1] > maxlat && bufy[n-2] > maxlat)) {
00157               n--;
00158            }
00159            else {
00160               break;
00161            }
00162        }
00163        if (n <= 1) continue;
00164 */
00165 
00166        if (mapform != NULL) (*mapform)(n,bufx,bufy); /* moved transformation to here   */
00167                                                      /* so bound-checking worked right */
00168 
00169        wrap = 0;
00170        /* check for wrap around problems */
00171        for (i = 0; i < n-1; i++) {
00172 
00173          /* jc: this code is wrong, as the bufx/y are floats that are
00174             converted to ints before abs() is called. Thus, small
00175             differences are masked off. The proof that it is wrong is
00176             that when replacing abs() with fabs(), as it should be,
00177             the code works the wrong way. What a proof :-)
00178 
00179            test[i] = abs((bufx[i]-bufx[i+1])) > abs(bufy[i]/3);
00180 
00181            If the intended behaviour is *really* that, than an
00182            explicit cast should be used to help other programmers, as
00183            is done bellow!!!
00184          */
00185 
00186            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" */
00187            if (test[i]) wrap = 1;
00188        }
00189 
00190        if (wrap == 0) {     
00191            plline(n,bufx,bufy);
00192        }
00193        else {
00194            for (i = 0; i < n-1; i++) {
00195               x[0] = bufx[i];
00196               x[1] = bufx[i+1];
00197               y[0] = bufy[i];
00198               y[1] = bufy[i+1];
00199               if (test[i] == 0) {
00200                   plline(2,x,y);
00201               }
00202               else {  /* this code seems to supercede the block commented out above */
00203                   /* segment goes off the edge */
00204                   sign = (x[1] > x[0]) ? 1 : -1;
00205                   x[1] -= sign * 360.0;
00206                   plline(2,x,y);
00207                   x[0] = bufx[i];
00208                   x[1] = bufx[i+1];
00209                   y[0] = bufy[i];
00210                   y[1] = bufy[i+1];
00211                   x[0] += sign * 360.0;
00212                   plline(2,x,y);
00213               }
00214            }
00215        }
00216     }
00217 }
00218 
00219 /*----------------------------------------------------------------------*\
00220  * void plmeridians(void (*mapform)(PLINT, PLFLT *, PLFLT *), 
00221  *                PLFLT dlong, PLFLT dlat, PLFLT minlong, PLFLT maxlong, 
00222  *                PLFLT minlat, PLFLT maxlat);
00223  *
00224  * Plot the latitudes and longitudes on the background.  The lines 
00225  * are plotted in the current color and line style.
00226  * 
00227  * mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
00228  * coordinate longitudes and latitudes to a plot coordinate system.  By
00229  * using this transform, we can change from a longitude, latitude
00230  * coordinate to a polar stereographic project, for example.  Initially,
00231  * x[0]..x[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
00232  * latitudes.  After the call to mapform(), x[] and y[] should be replaced
00233  * by the corresponding plot coordinates.  If no transform is desired,
00234  * mapform can be replaced by NULL.
00235  * 
00236  * dlat, dlong are the interval in degrees that the latitude and longitude
00237  * lines are to be plotted. 
00238  * 
00239  * minlong, maxlong are the values of the longitude on the left and right
00240  * side of the plot, respectively. The value of minlong must be less than
00241  * the values of maxlong, and the values of maxlong-minlong must be less
00242  * or equal to 360.
00243  * 
00244  * minlat, maxlat are the minimum and maximum latitudes to be plotted on
00245  * the background.  One can always use -90.0 and 90.0 as the boundary
00246  * outside the plot window will be automatically eliminated.  However, the
00247  * program will be faster if one can reduce the size of the background
00248  * plotted.
00249 \*----------------------------------------------------------------------*/
00250 
00251 #define NSEG 100
00252 
00253 void 
00254 plmeridians( void (*mapform)(PLINT, PLFLT *, PLFLT *), 
00255                PLFLT dlong, PLFLT dlat,
00256                PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
00257 {
00258     PLFLT yy, xx, temp, x[2], y[2], dx, dy;
00259 
00260     if (minlong > maxlong) {
00261        temp = minlong;
00262        minlong = maxlong;
00263        maxlong = temp;
00264     }
00265     if (minlat > maxlat) {
00266        temp = minlat;
00267        minlat = maxlat;
00268        maxlat = temp;
00269     }
00270     dx = (maxlong - minlong) / NSEG;
00271     dy = (maxlat - minlat) / NSEG;
00272 
00273     /* latitudes */
00274 
00275     for (yy = dlat * ceil(minlat/dlat); yy <= maxlat; yy += dlat) {
00276        if (mapform == NULL) {
00277            y[0] = y[1] = yy;
00278            x[0] = minlong;
00279            x[1] = maxlong;
00280            plline(2,x,y);
00281        }
00282        else {
00283            for (xx = minlong; xx < maxlong; xx += dx) {
00284                y[0] = y[1] = yy;
00285               x[0] = xx;
00286               x[1] = xx + dx;
00287               (*mapform)(2,x,y);
00288               plline(2,x,y);
00289            }
00290        }
00291     }
00292 
00293     /* longitudes */
00294  
00295     for (xx = dlong * ceil(minlong/dlong); xx <= maxlong; xx += dlong) {
00296         if (mapform == NULL) {
00297             x[0] = x[1] = xx;
00298             y[0] = minlat;
00299             y[1] = maxlat;
00300             plline(2,x,y);
00301         }
00302         else {
00303             for (yy = minlat; yy < maxlat; yy += dy) {
00304                 x[0] = x[1] = xx;
00305                 y[0] = yy;
00306                 y[1] = yy + dy;
00307                 (*mapform)(2,x,y);
00308                 plline(2,x,y);
00309             }
00310         }
00311     }
00312 }