plmap.c

Go to the documentation of this file.

//      Continental Outline and Political Boundary Backgrounds
//
//      Some plots need a geographical background such as the global
//      surface temperatures or the population density.  The routine
//      plmap() will draw one of the following backgrounds: continental
//      outlines, political boundaries, the United States, and the United
//      States with the continental outlines.  The routine plmeridians()
//      will add the latitudes and longitudes to the background.  After
//      the background has been drawn, one can use a contour routine or a
//      symbol plotter to finish off the plot.
//
//      Copyright (C) 1991, 1993, 1994  Wesley Ebisuzaki
//      Copyright (C) 1994, 2000, 2001  Maurice LeBrun
//      Copyright (C) 1999  Geoffrey Furnish
//      Copyright (C) 2000-2016 Alan W. Irwin
//      Copyright (C) 2001  Andrew Roach
//      Copyright (C) 2001, 2004  Rafael Laboissiere
//      Copyright (C) 2002  Vincent Darley
//      Copyright (C) 2003  Joao Cardoso
//
//      This file is part of PLplot.
//
//      PLplot is free software; you can redistribute it and/or modify
//      it under the terms of the GNU Library General Public License
//      as published by the Free Software Foundation; version 2 of the
//      License.
//
//      PLplot is distributed in the hope that it will be useful, but
//      WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU Library General Public License for more details.
//
//      You should have received a copy of the GNU Library General Public
//      License along with this library; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301
//      USA
//
//

#define DEBUG
#define NEED_PLDEBUG

#include "plplotP.h"

#ifdef HAVE_SHAPELIB
#include <shapefil.h>

SHPHandle
OpenShapeFile( PLCHAR_VECTOR fn );

#ifdef HAVE_SAHOOKS
static void
CustomErrors( PLCHAR_VECTOR message );
#endif //HAVE_SAHOOKS

#define OpenMap     OpenShapeFile
#define CloseMap    SHPClose

//redistributes the lon value so that it is within +/-180 deg of midlon
//for wrapping purposes.
void
rebaselon( PLFLT *lon, PLFLT midlon )
{
    if ( *lon > midlon + 180.0 )
        *lon -= floor( ( *lon - midlon - 180.0 ) / 360.0 + 1.0 ) * 360.0;
    else if ( *lon < midlon - 180.0 )
        *lon += floor( ( midlon - 180.0 - *lon ) / 360.0 + 1.0 ) * 360.0;
}

//--------------------------------------------------------------------------
//Actually draw the map lines points and text.
//--------------------------------------------------------------------------
void
drawmapdata( PLMAPFORM_callback mapform, int shapetype, PLINT n, double *x, double *y, PLFLT dx, PLFLT dy, PLFLT just, PLCHAR_VECTOR text )
{
    PLINT i;
    PLFLT *renderX;
    PLFLT *renderY;

    //we need to do something a bit different with filled polygons. The issue is the poles
    //on lat/lon plots. If we draw Antarctica then we expect it to be filled, but this
    //requires us to add in the extra points extending to the pole.
    //We identify a fill that loops round the globe because its start point will not be its
    //end point due to the processing done in drawmap. It should be 360 degrees different to
    //within the rounding error.
    //For a basic plot without any transformation we could just add three points to go down
    //to the pole accross to the starting longitude then up to join the start point, but if
    //we consider a polar plot looking down on the North Pole, then Antarctica gets turned
    //inside out so we actually want to add many points at the pole with different longitudes.

    if ( ( shapetype == SHPT_POLYGON || shapetype == SHPT_POLYGONZ || shapetype == SHPT_POLYGONM ) && x[0] != x[n - 1] )
    {
        //we have a polygon that rings the pole - render it differently to everything else

        if ( x[0] != x[n - 1] )
        {
            //this is a looped round the pole polygon
            PLINT  i;
            double poleLat;
            double lonInterval;
            PLINT  nExtraPoints = MAX( 501, n + 1 );
            double *newX;
            double *newY;

            //The shapefile standard says that as we traverse the vertices, the inside should be
            //on the right
            if ( x[n - 1] > x[0] )
                poleLat = -90.0;
            else
                poleLat = 90.0;

            newX = malloc( ( n + nExtraPoints ) * sizeof ( double ) );
            newY = malloc( ( n + nExtraPoints ) * sizeof ( double ) );
            if ( !newX || !newY )
            {
                free( newX );
                free( newY );
                plabort( "Could not allocate memory for adding pole points to a map polygon." );
                return;
            }
            memcpy( newX, x, n * sizeof ( double ) );
            memcpy( newY, y, n * sizeof ( double ) );

            lonInterval = ( x[0] - x[n - 1] ) / (double) ( nExtraPoints - 2 );

            for ( i = 0; i < nExtraPoints - 1; ++i )
            {
                newX[n + i] = x[n - 1] + i * lonInterval;
                newY[n + i] = poleLat;
            }
            newX[n + nExtraPoints - 1] = x[0];
            newY[n + nExtraPoints - 1] = y[0];
#ifdef PL_DOUBLE
            renderX = newX;
            renderY = newY;
#else
            fltX = malloc( ( n + nExtraPoints ) * sizeof ( PLFLT ) );
            fltX = malloc( ( n + nExtraPoints ) * sizeof ( PLFLT ) );
            if ( !fltX || !fltY )
            {
                free( fltX );
                free( fltY );
                free( newX );
                free( newY );
                plabort( "Could not allocate memory converting map date to floats." );
                return;
            }
            for ( i = 0; i < n + nExtraPoints; ++i )
            {
                fltX[i] = (PLFLT) newX[i];
                fltY[i] = (PLFLT) newY[i];
            }
            renderX = fltX;
            renderY = fltY;
#endif
            if ( mapform != NULL )
                (*mapform)( n + nExtraPoints, renderX, renderY );
            plfill( n + nExtraPoints, renderX, renderY );
            free( newX );
            free( newY );
#ifndef PL_DOUBLE
            free( fltX );
            free( fltY );
#endif

            //rendering is complete, return;
            return;
        }
        else
        {
            //this is just a regular polygon - render it as we would expect
            if ( mapform != NULL )
                (*mapform)( n, x, y );
            plfill( n, x, y );
        }

        //return here as we are done
        return;
    }

    //if we get to here we don't have a polygon wrapping all the way round the globe
    //Just do normal rendering

    //convert to floats if needed
#ifdef PL_DOUBLE
    renderX = x;
    renderY = y;
#else
    fltX = malloc( n * sizeof ( PLFLT ) );
    fltX = malloc( n * sizeof ( PLFLT ) );
    if ( !fltX || !fltY )
    {
        free( fltX );
        free( fltY );
        plabort( "Could not allocate memory converting map date to floats." );
        return;
    }
    for ( i = 0; i < n; ++i )
    {
        fltX[i] = (PLFLT) x[i];
        fltY[i] = (PLFLT) y[i];
    }
    renderX = fltX;
    renderY = fltY;
#endif

    //do the transform if needed
    if ( mapform != NULL )
        ( *mapform )( n, renderX, renderY );

    //deo the rendering
    if ( shapetype == SHPT_ARC )
        plline( n, renderX, renderY );
    else if ( shapetype == SHPT_POINT )
        for ( i = 0; i < n; ++i )
            plptex( renderX[i], renderY[i], dx, dy, just, text );
    else if ( shapetype == SHPT_ARCZ || shapetype == SHPT_ARCM )
        plline( n, renderX, renderY );
    else if ( shapetype == SHPT_POINTM || shapetype == SHPT_POINTZ )
        for ( i = 0; i < n; ++i )
            plptex( renderX[i], renderY[i], dx, dy, just, text );
    else if ( shapetype == SHPT_POLYGON || shapetype == SHPT_POLYGONZ || shapetype == SHPT_POLYGONM )
        plfill( n, renderX, renderY );
    else if ( shapetype == SHPT_NULL )
        ;         //do nothing for a NULL shape
    else
        plabort( "Unknown render type passed in to drawmapdata(). PLplot can render shapefile arcs, points and polygons (including z and m variants)." );

#ifndef PL_DOUBLE
    free( fltX );
    free( fltY );
#endif
}


//--------------------------------------------------------------------------
//This is a function called by the front end map functions to do the map drawing. Its
//parameters are:
//mapform: The transform used to convert the data in raw coordinates to x, y positions
//on the plot
//name: either one of the plplot provided lat/lon maps or the path/file name of a
//shapefile
//dx/dy: the gradient of text/symbols drawn if text is non-null
//shapetype: one of ARC, SHPT_ARCZ, SHPT_ARCM, SHPT_POLYGON, SHPT_POLYGONZ,
//SHPT_POLYGONM, SHPT_POINT, SHPT_POINTM, SHPT_POINTZ. See drawmapdata() for the
//how each type is rendered. But Basically the ARC options are lines, the POLYGON
//options are filled polygons, the POINT options are points/text. Options beginning
//SHPT will only be defined if HAVE_SHAPELIB is true
//text: The text (which can be actual text or a unicode symbol) to be drawn at
//each point
//minx/maxx: The min/max longitude when using a plplot provided map or x value if
//using a shapefile
//miny/maxy: The min/max latitude when using a plplot provided map or y value if
//using a shapefile
//plotentries: used only for shapefiles, as one shapefile contains multiple vectors
//each representing a different item (e.g. multiple boundaries, multiple height
//contours etc. plotentries is an array containing the indices of the
//entries within the shapefile that you wish to plot. if plotentries is null all
//entries are plotted
//nplotentries: the number of elements in plotentries. Ignored if plplot was not built
//with shapefile support or if plotentries is null
//--------------------------------------------------------------------------
void
drawmap( PLMAPFORM_callback mapform, PLCHAR_VECTOR name,
         PLFLT dx, PLFLT dy, int shapetype, PLFLT just, PLCHAR_VECTOR text,
         PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy, PLINT_VECTOR plotentries, PLINT nplotentries )
{
    int    i;
    char   *filename = NULL;                               //filename of the map data
    char   warning[1024];                                  //string used for constructing a sensible warning message
    int    nVertices = 200;                                //number of vertices in a particular part
    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
    //PLFLT  *bufx   = NULL, *bufy = NULL; //the data that will be plot after dealing with wrapping round the globe
    int    bufsize = 0;                                    //number of elements in bufx and bufy
    size_t filenamelen;                                    //length of the filename
    char   islatlon     = 1;                               //specify if the data is lat-lon or some other projection
    int    appendresult = 0;


    SHPHandle in;
    int       nentries;                            //number of objects in the shapefile
    int       entryindex  = 0;                     //index of plotentries that we are currently rendering
    int       entrynumber = 0;                     //id of the object we are currently rendering
    int       partnumber  = 0;                     //part of the object we are currently rendering (some objects are split into parts)
    SHPObject *object     = NULL;                  //pointer to the object we are currently rendering
    double    *bufxraw;                            //pointer to the raw x data read from the file
    double    *bufyraw;                            //pointer to the raw y data read from the file
    char      *prjfilename = NULL;                 //filename of the projection file
    PDFstrm   *prjfile;                            //projection file
    char      prjtype[] = { 0, 0, 0, 0, 0, 0, 0 }; //string holding the projection type description
    double    fileMins[4];                         //min x, y, z, m for the file
    double    fileMaxs[4];                         //max x, y, z, m for the file
    int       fileShapeType;                       //the shapetype read from the file

    //
    // read map outline
    //

    //strip the .shp extension if a shapefile has been provided
    if ( strstr( name, ".shp" ) )
        filenamelen = ( strlen( name ) - 4 );
    else
        filenamelen = strlen( name );
    filename = (char *) malloc( filenamelen + 1 );
    if ( !filename )
    {
        plabort( "Could not allocate memory for map filename root" );
        return;
    }
    strncpy( filename, name, filenamelen );
    filename[ filenamelen ] = '\0';

    //Open the shp and shx file using shapelib
    if ( ( in = OpenShapeFile( filename ) ) == NULL )
    {
        strcpy( warning, "Could not find " );
        strcat( warning, filename );
        strcat( warning, " file." );
        plabort( warning );
        free( filename );
        return;
    }
    SHPGetInfo( in, &nentries, &fileShapeType, fileMins, fileMaxs );
    if ( shapetype == SHPT_NULL )
    {
        shapetype = fileShapeType;
    }
    //also check for a prj file which will tell us if the data is lat/lon or projected
    //if it is projected then set ncopies to 1 - i.e. don't wrap round longitudes
    prjfilename = (char *) malloc( filenamelen + 5 );
    if ( !prjfilename )
    {
        free( filename );
        plabort( "Could not allocate memory for generating map projection filename" );
        return;
    }
    strncpy( prjfilename, name, filenamelen );
    prjfilename[ filenamelen ] = '\0';
    strcat( prjfilename, ".prj" );
    prjfile = plLibOpenPdfstrm( prjfilename );
    if ( prjfile && prjfile->file )
    {
        fread( prjtype, 1, 6, prjfile->file );
        if ( strcmp( prjtype, "PROJCS" ) == 0 )
            islatlon = 0;
        pdf_close( prjfile );
    }
    free( prjfilename );
    prjfilename = NULL;

    for (;; )
    {
        //each object in the shapefile is split into parts.
        //If we are need to plot the first part of an object then read in a new object
        //and check how many parts it has. Otherwise use the object->panPartStart vector
        //to check the offset of this part and the next part and allocate memory. Copy
        //the data to this memory converting it to PLFLT and draw it.
        //finally increment the part number or if we have finished with the object reset the
        //part numberand increment the object.

        //break condition if we've reached the end of the file
        if ( ( !plotentries && ( entrynumber == nentries ) ) || ( plotentries && ( entryindex == nplotentries ) ) )
            break;

        //if partnumber == 0 then we need to load the next object
        if ( partnumber == 0 )
        {
            if ( plotentries )
                object = SHPReadObject( in, plotentries[entryindex] );
            else
                object = SHPReadObject( in, entrynumber );
        }
        //if the object could not be read, increment the object index to read and
        //return to the top of the loop to try the next object.
        if ( object == NULL )
        {
            entrynumber++;
            entryindex++;
            partnumber = 0;
            continue;
        }

        //work out how many points are in the current part
        if ( object->nParts == 0 )
            nVertices = object->nVertices;                                                       //if object->nParts==0, we can still have 1 vertex. A bit odd but it's the way it goes
        else if ( partnumber == ( object->nParts - 1 ) )
            nVertices = object->nVertices - object->panPartStart[partnumber];                    //panPartStart holds the offset for each part
        else
            nVertices = object->panPartStart[partnumber + 1] - object->panPartStart[partnumber]; //panPartStart holds the offset for each part

        //point the plot buffer to the correct starting vertex
        //and copy it to the PLFLT arrays. If we had object->nParts == 0
        //then panPartStart will be NULL
        if ( object->nParts > 0 )
        {
            bufxraw = object->padfX + object->panPartStart[partnumber];
            bufyraw = object->padfY + object->panPartStart[partnumber];
        }
        else
        {
            bufxraw = object->padfX;
            bufyraw = object->padfY;
        }

        //set the min/max x/y of the object
        minsectlon = object->dfXMin;
        maxsectlon = object->dfXMax;
        minsectlat = object->dfYMin;
        maxsectlat = object->dfYMax;

        if ( nVertices > 0 )
        {
            if ( islatlon )
            {
                PLINT  nExtraPositiveRedraws = 0;
                double unrebasedFirstValue   = bufxraw[0];
                double rebaseAmount;
                //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
                //   or even more bizzare options - like google earth when zoomed
                //   right out loops round and round.
                //
                //we can cure these problems by conditionally adding/subtracting
                //360 degrees to each longitude point in order to ensure that the
                //distance between adgacent points is always less than 180
                //degrees, then repeatedly shifting the shape by integer multiples
                //of 360 in each longitude direction until the shape has moved out
                //of the plot area.

                //reset our plot longitude limits to avoid crazy loops in case the user has
                //set them to be infinities or similar
                if ( minx == -PLFLT_MAX || minx <= -PLFLT_HUGE_VAL )
                    minx = fileMins[0];
                if ( maxx == -PLFLT_MAX || maxx >= PLFLT_HUGE_VAL )
                    maxx = fileMaxs[0];

                //move the first point by multiples of 360 putting it as close
                //as possible to our plot limits centre point.
                unrebasedFirstValue = bufxraw[0];
                rebaselon( &bufxraw[0], ( minx + maxx ) / 2.0 );
                rebaseAmount = bufxraw[0] - unrebasedFirstValue;
                minsectlon  += rebaseAmount;
                maxsectlon  += rebaseAmount;
                for ( i = 1; i < nVertices; i++ )
                {
                    double difference;
                    //first adjust the point by the same amount as our first point
                    bufxraw[i] += rebaseAmount;

                    //now check it doesn't do any strange wrapping
                    difference = bufxraw[i] - bufxraw[i - 1];
                    if ( difference < -180.0 )
                    {
                        bufxraw[i] += 360.0;
                        maxsectlon  = MAX( maxsectlon, bufxraw[i] );
                    }
                    else if ( difference > 180.0 )
                    {
                        bufxraw[i] -= 360.0;
                        minsectlon  = MIN( minsectlon, bufxraw[i] );
                    }
                }

                //check if the latitude and longitude range means we need to plot this section
                if ( ( maxsectlat > miny ) && ( minsectlat < maxy )
                     && ( maxsectlon > minx ) && ( minsectlon < maxx ) )
                {
                    drawmapdata( mapform, shapetype, nVertices, bufxraw, bufyraw, dx, dy, just, text );
                }
                //now check if adding multiples of 360 to the data still allow it to be plotted
                while ( minsectlon + 360.0 < maxx )
                {
                    for ( i = 0; i < nVertices; ++i )
                        bufxraw[i] += 360.0;
                    drawmapdata( mapform, shapetype, nVertices, bufxraw, bufyraw, dx, dy, just, text );
                    minsectlon += 360.0;
                    ++nExtraPositiveRedraws;
                }
                if ( maxsectlon - 360.0 > minx )
                {
                    for ( i = 0; i < nVertices; ++i )
                        bufxraw[i] -= nExtraPositiveRedraws * 360.0;
                    while ( maxsectlon - 360.0 > minx )
                    {
                        for ( i = 0; i < nVertices; ++i )
                            bufxraw[i] -= 360.0;
                        drawmapdata( mapform, shapetype, nVertices, bufxraw, bufyraw, dx, dy, just, text );
                        maxsectlon -= 360.0;
                    }
                }
            }
            else
            {
                drawmapdata( mapform, shapetype, nVertices, bufxraw, bufyraw, dx, dy, just, text );
            }
        }


        //increment the partnumber or if we've reached the end of
        //an entry increment the entrynumber and set partnumber to 0
        if ( partnumber == object->nParts - 1 || object->nParts == 0 )
        {
            entrynumber++;
            entryindex++;
            partnumber = 0;
            SHPDestroyObject( object );
            object = NULL;
        }
        else
            partnumber++;
    }
    // Close map file
    SHPClose( in );

    //free memory
    free( filename );
}
#endif //HAVE_SHAPELIB


//--------------------------------------------------------------------------
// void plmap(PLMAPFORM_callback mapform, PLCHAR_VECTOR name,
//            PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy);
//
// plot continental outline in world coordinates
//
// v1.4: machine independant version
// v1.3: replaced plcontinent by plmap, added plmeridians
// v1.2: 2 arguments:  mapform, type of plot
//
// mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
// coordinate longitudes and latitudes to a plot coordinate system.  By
// using this transform, we can change from a longitude, latitude
// coordinate to a polar stereographic project, for example.  Initially,
// x[0]..[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
// latitudes.  After the call to mapform(), x[] and y[] should be replaced
// by the corresponding plot coordinates.  If no transform is desired,
// mapform can be replaced by NULL.
//
// name is a character string. The value of this parameter determines the
// name of background. The possible values are,
//
//      "globe"         continental outlines
//      "usa"           USA and state boundaries
//      "cglobe"        continental outlines and countries
//      "usaglobe"      USA, state boundaries and continental outlines
// alternatively the filename of a shapefile can be passed if PLplot has
// been compiled with shapelib. In this case either the base name of the
// file can be passed or the filename including the .shp or .shx suffix.
// Only the .shp and .shx files are used.
//
// minx, maxx are the minimum and maximum untransformed x values to be
// plotted. For the built in plots these are longitudes. For other
//shapefiles these are in the units of the shapefile. The value of minx
//must be less than the values of maxx.
//
// miny, maxy are as per minx and maxx except for the built in plots
//the units are latitude.
//--------------------------------------------------------------------------

void
plmap( PLMAPFORM_callback mapform, PLCHAR_VECTOR name,
       PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy )
{
#ifdef HAVE_SHAPELIB
    drawmap( mapform, name, 0.0, 0.0, SHPT_NULL, 0.0, NULL, minx, maxx,
        miny, maxy, NULL, 0 );
#else
    plwarn( "plmap is a no-op because shapelib is not available." );
#endif
}

//--------------------------------------------------------------------------
// void plmapline(PLMAPFORM_callback mapform,
//              PLCHAR_VECTOR name, PLFLT minx, PLFLT maxx, PLFLT miny,
//              PLFLT maxy, PLINT_VECTOR plotentries, PLINT nplotentries);

//New version of plmap which allows us to specify which items in a shapefile
//we want to use. parameters are as above but with the plotentries being an
//array containing the indices of the elements we wish to draw and
//nplotentries being the number of items in plotentries.
//If shapefile access was not built into plplot then plotentries and
//nplotentries are ignored. If plotentries is null than all entries are
//drawn and nplotentries is ignored.
//The name distiguishes it from other functions which plot points/text and
//polygons, but note that the type of element in the shapefile need not
//match the type of element drawn - i.e. arc elements from a shapefile could
//be drawn as points using the plmaptex function.
//--------------------------------------------------------------------------
void
plmapline( PLMAPFORM_callback mapform, PLCHAR_VECTOR name,
           PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy,
           PLINT_VECTOR plotentries, PLINT nplotentries )
{
#ifdef HAVE_SHAPELIB
    drawmap( mapform, name, 0.0, 0.0, SHPT_ARC, 0.0, "", minx, maxx,
        miny, maxy, plotentries, nplotentries );
#else
    plwarn( "plmapline is a no-op because shapelib is not available." );
#endif
}

//--------------------------------------------------------------------------
// void plmapstring(PLMAPFORM_callback mapform,
//              PLCHAR_VECTOR name, PLFLT just, PLCHAR_VECTOR string,
//              PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy,
//              PLINT_VECTOR plotentries, PLINT nplotentries);
//
//As per plmapline but plots symbols. The map equivalent of plstring. string
//has the same meaning as in plstring.
//--------------------------------------------------------------------------
void
plmapstring( PLMAPFORM_callback mapform,
             PLCHAR_VECTOR name, PLCHAR_VECTOR string,
             PLFLT minx, PLFLT maxx, PLFLT miny, PLFLT maxy,
             PLINT_VECTOR plotentries, PLINT nplotentries )
{
#ifdef HAVE_SHAPELIB
    drawmap( mapform, name, 1.0, 0.0, SHPT_POINT, 0.5, string, minx, maxx,
        miny, maxy, plotentries, nplotentries );
#else
    plwarn( "plmapstring is a no-op because shapelib is not available." );
#endif
}

//--------------------------------------------------------------------------
// 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);
//
//As per plmapline but plots text. The map equivalent of plptex. dx, dy,
//just and text have the same meaning as in plptex.
//--------------------------------------------------------------------------
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 )
{
#ifdef HAVE_SHAPELIB
    drawmap( mapform, name, dx, dy, SHPT_POINT, just, text, minx, maxx,
        miny, maxy, &plotentry, 1 );
#else
    plwarn( "plmaptex is a no-op because shapelib is not available." );
#endif
}

//--------------------------------------------------------------------------
// void plmapfill(PLMAPFORM_callback mapform,
//              PLCHAR_VECTOR name, PLFLT minx, PLFLT maxx, PLFLT miny,
//              PLFLT maxy, PLINT_VECTOR plotentries, PLINT nplotentries);
//
//As per plmapline but plots a filled polygon. The map equivalent to
//plfill. Uses the pattern defined by plsty or plpat.
//--------------------------------------------------------------------------
void
plmapfill( PLMAPFORM_callback mapform,
           PLCHAR_VECTOR name, PLFLT minx, PLFLT maxx, PLFLT miny,
           PLFLT maxy, PLINT_VECTOR plotentries, PLINT nplotentries )
{
#ifdef HAVE_SHAPELIB
    drawmap( mapform, name, 0.0, 0.0, SHPT_POLYGON, 0.0, NULL, minx, maxx,
        miny, maxy, plotentries, nplotentries );
#else
    plwarn( "plmapfill is a no-op because shapelib is not available." );
#endif
}

//--------------------------------------------------------------------------
// void plmeridians(PLMAPFORM_callback mapform,
//                  PLFLT dlong, PLFLT dlat, PLFLT minx, PLFLT maxx,
//                  PLFLT miny, PLFLT maxy);
//
// Plot the latitudes and longitudes on the background.  The lines
// are plotted in the current color and line style.
//
// mapform(PLINT n, PLFLT *x, PLFLT *y) is a routine to transform the
// coordinate longitudes and latitudes to a plot coordinate system.  By
// using this transform, we can change from a longitude, latitude
// coordinate to a polar stereographic project, for example.  Initially,
// x[0]..x[n-1] are the longitudes and y[0]..y[n-1] are the corresponding
// latitudes.  After the call to mapform(), x[] and y[] should be replaced
// by the corresponding plot coordinates.  If no transform is desired,
// mapform can be replaced by NULL.
//
// dlat, dlong are the interval in degrees that the latitude and longitude
// lines are to be plotted.
//
// minx, maxx are the values of the longitude on the left and right
// side of the plot, respectively. The value of minx must be less than
// the values of maxx, and the values of maxx-minx must be less
// or equal to 360.
//
// miny, maxy are the minimum and maximum latitudes to be plotted on
// the background.  One can always use -90.0 and 90.0 as the boundary
// outside the plot window will be automatically eliminated.  However, the
// program will be faster if one can reduce the size of the background
// plotted.
//--------------------------------------------------------------------------

#define NSEG    100

void
plmeridians( PLMAPFORM_callback mapform,
             PLFLT dlong, PLFLT dlat,
             PLFLT minlong, PLFLT maxlong, PLFLT minlat, PLFLT maxlat )
{
    PLFLT yy, xx, temp, x[2], y[2], dx, dy;

    if ( minlong > maxlong )
    {
        temp    = minlong;
        minlong = maxlong;
        maxlong = temp;
    }
    if ( minlat > maxlat )
    {
        temp   = minlat;
        minlat = maxlat;
        maxlat = temp;
    }
    dx = ( maxlong - minlong ) / NSEG;
    dy = ( maxlat - minlat ) / NSEG;

    // latitudes

    for ( yy = dlat * ceil( minlat / dlat ); yy <= maxlat; yy += dlat )
    {
        if ( mapform == NULL )
        {
            plpath( NSEG, minlong, yy, maxlong, yy );
        }
        else
        {
            for ( xx = minlong; xx < maxlong; xx += dx )
            {
                y[0] = y[1] = yy;
                x[0] = xx;
                x[1] = xx + dx;
                ( *mapform )( 2, x, y );
                plline( 2, x, y );
            }
        }
    }

    // longitudes

    for ( xx = dlong * ceil( minlong / dlong ); xx <= maxlong; xx += dlong )
    {
        if ( mapform == NULL )
        {
            plpath( NSEG, xx, minlat, xx, maxlat );
        }
        else
        {
            for ( yy = minlat; yy < maxlat; yy += dy )
            {
                x[0] = x[1] = xx;
                y[0] = yy;
                y[1] = yy + dy;
                ( *mapform )( 2, x, y );
                plline( 2, x, y );
            }
        }
    }
}

//--------------------------------------------------------------------------
// SHPHandle OpenShapeFile(fn)
//
//--------------------------------------------------------------------------
#ifdef HAVE_SHAPELIB
#ifdef HAVE_SAHOOKS
// Our thanks to Frank Warmerdam, the developer of shapelib for suggesting
// this approach for quieting shapelib "Unable to open" error messages.
static
void CustomErrors( PLCHAR_VECTOR message )
{
    if ( strstr( message, "Unable to open" ) == NULL )
        fprintf( stderr, "%s\n", message );
}
#endif

SHPHandle
OpenShapeFile( PLCHAR_VECTOR fn )
{
    SHPHandle file;
    char      *fs = NULL, *dn = NULL;
#ifdef HAVE_SAHOOKS
    SAHooks   sHooks;

    SASetupDefaultHooks( &sHooks );
    sHooks.Error = CustomErrors;
#else
    // Using ancient version of shapelib without SAHooks or SHPOpenLL.
    // For this case live with the misleading "Unable to open" error
    // messages.
    // int sHooks;
#define SHPOpenLL( a, b, c )    SHPOpen( a, b )
#endif //HAVE_SAHOOKS

//***   search build tree               ***

    if ( plInBuildTree() == 1 )
    {
        plGetName( SOURCE_DIR, "data", fn, &fs );

        if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != 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 = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
            goto done;
        fprintf( stderr, PLPLOT_LIB_ENV "=\"%s\"\n", dn ); // what IS set?
    }
#endif  // PLPLOT_LIB_ENV

//***   search current directory        ***

    if ( ( file = SHPOpenLL( fn, "rb", &sHooks ) ) != NULL )
    {
        pldebug( "OpenShapeFile", "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 = SHPOpenLL( fs, "rb", &sHooks ) ) != 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 = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
        goto done;
#endif  // DATA_DIR

//***   search hardwired location       ***

#ifdef PLLIBDEV
    plGetName( PLLIBDEV, "", fn, &fs );

    if ( ( file = SHPOpenLL( fs, "rb", &sHooks ) ) != NULL )
        goto done;
#endif  // PLLIBDEV

//***   not found, give up      ***
    pldebug( "OpenShapeFile", "File %s not found.\n", fn );
    free_mem( fs );
    return NULL;

done:
    pldebug( "OpenShapeFile", "SHPOpen successfully opened two files with basename %s\n", fs );
    free_mem( fs );
    return ( file );
}
#endif //ifdef HAVE_SHAPELIB