Logo Search packages:      
Sourcecode: magics++ version File versions  Download package

gshhs2magics.cc

/*    $Id: gshhs_dp.c,v 1.5 2004/09/28 18:09:25 pwessel Exp $
 *
 * gshhs_dp applies the Douglas-Peucker algorithm to simplify a line
 * segment given a tolerance.  The algorithm is based on the paper
 * Douglas, D. H., and T. K. Peucker, Algorithms for the reduction
 *   of the number of points required to represent a digitized line
 *   of its caricature, Can. Cartogr., 10, 112-122, 1973.
 * The implementation of this algorithm has been kindly provided by
 * Dr. Gary J. Robinson, Environmental Systems Science Centre,
 * University of Reading, Reading, UK (gazza@mail.nerc-essc.ac.uk)
 *
 * Paul Wessel, 18-MAY-1999
 * Version: 1.1 Added byte flipping
 *        1.2 Explicit binary read for DOS.  POSIX compliance
 *        1.3, 08-NOV-1999: Released under GNU GPL
 *        1.4 05-SEPT-2000: Made a GMT supplement; FLIP no longer needed
 *        1.5 11-SEPT-2004: Updated to work with GSHHS v1.3 data format
 *
 *    Copyright (c) 1996-2004 by P. Wessel and W. H. F. Smith
 *    See COPYING file for copying and redistribution conditions.
 *
 *    This program is free software; you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation; version 2 of the License.
 *
 *    This program 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 General Public License for more details.
 *
 *    Contact info: www.soest.hawaii.edu/pwessel
 */

using namespace std;

#include <fstream>

#include "gshhs.h"
#include "Polyline.h"
#include "GeoPoint.h"
#include "netcdfcpp.h"
#include "Netcdf.h"


#define sqr(x) ((x)*(x))
#define D2R (M_PI/180.0)
#define F (D2R * 0.5 * 1.0e-6)
#define FALSE 0
#define VNULL (void *)NULL

int Douglas_Peucker_i (int x_source[], int y_source[], int n_source, double band, int index[]);
void *get_memory (void *prev_addr, int n, size_t size, char *progname);


template <class P>
double area(Polyline& c)
{
            double area = 0 ; 
            if(c.size()) c.push_back(c[0]);

            for(unsigned int j = 0; j < c.size(); j++)
            { 
                  if(j) area += (c[j-1].x() * c[j].y() - c[j].x() * c[j-1].y() );
            }
            if(c.size()) c.pop_back();
            return area;
}

int main (int argc, char **argv)
{
      FILE  *fp_in;//, *fp_out;
      int   n_id, n_out, n=0, k, verbose = FALSE, *x, *y, *index;
      int   n_tot_in, n_tot_out, n_use, n_read, flip;
      double      redux, redux2, tolerance = 0.0;
      struct      GSHHS h;
      struct      POINT p;


      if (argc < 2 || !(argc == 6 || argc == 7)) {
            fprintf (stderr, "gshhs_magics v. 1.5 Line reduction using the Douglas-Peucker algorithm\n\n");
            fprintf (stderr, "usage:  gshhs_magics input.b tolerance <min land area Km**2> <min lake area Km**2>output.b [-v]\n");
            fprintf (stderr, "\ttolerance is maximum mismatch in km\n");
            fprintf (stderr, "\t-v will run in verbose mode and report shrinkage\n");
            exit (EXIT_FAILURE);
      }

      verbose = (argc == 7);
      tolerance = atof (argv[2]);
      if (verbose) fprintf (stderr,"gshhs_magics: Tolerance used is %g km\n", tolerance);
      fp_in  = fopen(argv[1], "rb");

      if ( !fp_in ) {
            fprintf (stderr,"gshhs_magics:  ERROR opening file %s.\n", argv[1]);
            exit(1);
      }

      string out = argv[5];
      
      /* Start shrink loop */
      
      n_id = n_out = n_tot_in = n_tot_out = 0;
      
      x = (int *) get_memory (VNULL, 1, sizeof (int), "gshhs_magics");
      y = (int *) get_memory (VNULL, 1, sizeof (int), "gshhs_magics");
      index = (int *) get_memory (VNULL, 1, sizeof (int), "gshhs_magics");
      
      n_read = fread ((void *)&h, sizeof (struct GSHHS), (size_t)1, fp_in);
      flip = (! (h.level > 0 && h.level < 5));

      vector<pair<struct GSHHS, Polyline* >  > coastlines;

      while (n_read == 1) {
      
            if (flip) {
                  h.id = swabi4 ((unsigned int)h.id);
                  h.n  = swabi4 ((unsigned int)h.n);
                  h.level = swabi4 ((unsigned int)h.level);
                  h.west  = swabi4 ((unsigned int)h.west);
                  h.east  = swabi4 ((unsigned int)h.east);
                  h.south = swabi4 ((unsigned int)h.south);
                  h.north = swabi4 ((unsigned int)h.north);
                  h.area  = swabi4 ((unsigned int)h.area);
                  h.version   = swabi4 ((unsigned int)h.version);
                  h.greenwich = swabi2 ((unsigned int)h.greenwich);
                  h.source    = swabi2 ((unsigned int)h.source);
        }

            if (verbose) fprintf (stderr, "Poly %6d", h.id);      
            
            x = (int *) get_memory ((void *)x, h.n, sizeof (int), "gshhs2magics");
            y = (int *) get_memory ((void *)y, h.n, sizeof (int), "gshhs2magics");
            index = (int *) get_memory ((void *)index, h.n, sizeof (int), "gshhs_magics");
            
            for (k = 0; k < h.n; k++) {
                  if (fread ((void *)&p, sizeof(struct POINT), (size_t)1, fp_in) != 1) {
                        fprintf (stderr,"gshhs_magics:  ERROR reading data point.\n");
                        exit (EXIT_FAILURE);
                  }
                  if (flip) {
                        p.x = swabi4 ((unsigned int)p.x);
                        p.y = swabi4 ((unsigned int)p.y);
                  }
                  x[k] = p.x;
                  y[k] = p.y;
            }
            int min_land_area = 10*atoi(argv[3]);
            int min_sea_area  = 10*atoi(argv[4]);
    
            
//        cout <<  "[" << level << "]";
//        cout <<  "[" << h.greenwich << "]";
//        cout << endl << h.area << "[" << level << "," << h.flag << "] ?? [" << min_land_area << ", " << min_sea_area << "]";

        if ( ((h.level == 1 ) && h.area >= min_land_area ) ||
              (h.level > 1    && h.area >= min_sea_area ) )
        {
                n_tot_in += h.n;

      //        n_use = (x[0] == x[h.n-1] && y[0] == y[h.n-1]) ? h.n-1 : h.n;
            n_use = h.n;

                //n = Douglas_Peucker_i (x, y, n_use, tolerance, index);

                if (n_use > 2) {
                      n++;
                      redux = 100.0 * (double) n / (double) h.n;
                      h.id = n_out;
                      h.n = n;
                Polyline*  poly = new Polyline();

                double lat, lon;

                for (k = 0; k < n_use; k++) {
                    lon = double(x[k])/1000000;
                    lat = double(y[k])/1000000;
                    //if ( k == 0 ) cout << lon << ", " << lat << endl;

                    poly->push_back(PaperPoint(lon, lat));
                      }


                if (tolerance > 0.0)
                {
                    if (n_use > 10)  // do not reduce polygons with few points
                        poly = poly->simplify(tolerance);
                    else if (n_use > 5)
                        poly = poly->simplify(tolerance/5.0);
                }


                coastlines.push_back(make_pair(h, poly));
                      n_out++;
                      n_tot_out += poly->size();
                }
                else
                      redux = 0.0;

                if (verbose) fprintf (stderr, "\t%.1f %% retained\n", redux);

                n_id++;
        }
        n_read = fread ((void *)&h, sizeof (struct GSHHS), (size_t)1, fp_in);
      }
   
    // Now we build the Netcdf file! 
    {
            NcFile netcdf = NcFile(out.c_str(), NcFile::Replace);
            NcDim* polys = netcdf.add_dim("polygones",coastlines.size());

            const long total_long = static_cast<long>(n_tot_out);

            NcDim* points = netcdf.add_dim("points", total_long );
            NcVar* position = netcdf.add_var("position", ncFloat, polys);
            NcVar* level = netcdf.add_var("level", ncFloat, polys);
            NcVar* area = netcdf.add_var("area", ncFloat, polys);
            NcVar* size = netcdf.add_var("size", ncFloat, polys);
            NcVar* greenwich = netcdf.add_var("greenwich", ncFloat, polys);
            NcVar* latitude = netcdf.add_var("latitude", ncFloat, points);
            NcVar* longitude = netcdf.add_var("longitude", ncFloat, points);
       
        
            int start = 0;
            int nb;
            int p = 0;
            vector<float> pos(coastlines.size());
            vector<float> nbs(coastlines.size());
            vector<float> levels(coastlines.size());
            vector<float> areas(coastlines.size());
            vector<float> greenwichs(coastlines.size());
            vector<float> lat(total_long);
            vector<float> lon(total_long);
            vector<int> pol(total_long);
            vector<int> p_v(total_long);
            int i = 0;
      
            for (vector<pair<struct GSHHS, Polyline* > >::const_iterator poly = coastlines.begin(); poly !=  coastlines.end(); ++poly)
            {
                  nb = poly->second->size();
                  pos[p] = start;
                  nbs[p] = nb;
            levels[p] = poly->first.level;
            areas[p] = poly->first.area;
            greenwichs[p] = poly->first.greenwich;
                  p++;
            //cout << "Area->" << poly->first.area << "(" << (*(poly->second)).area() << ") level->" << poly->first.level << "\n";
                  poly->second->makeClockwise();
                
            //cout << "revert now---> " << area(*(poly->second)) << "\n";
            
            for (Polyline::const_iterator point = poly->second->begin(); point !=  poly->second->end(); ++point) 
                  {
                        lat[i] =(*point).y();
                        lon[i] = (*point).x();
                //if (i ==0 ) 
                //    cout << "point" << lat[i] << ", " << lon[i] << endl;
                        i++;
                  }
                  start += nb; 
            }
       
            position->put(&pos[0], pos.size());
            size->put(&nbs[0], nbs.size());     
            level->put(&levels[0], levels.size());
            area->put(&areas[0], areas.size());
            greenwich->put(&greenwichs[0], greenwichs.size());
            latitude->put(&lat[0], lat.size());
            longitude->put(&lon[0], lon.size());
       
            netcdf.close();
        }
    
      free ((void *)x); 
      free ((void *)y); 
      free ((void *)index);   
            
      fclose (fp_in);

      redux = 100.0 * (1.0 - (double) n_tot_out / (double) n_tot_in);
      redux2 = 100.0 * (1.0 - (double) n_out / (double) n_id);
      printf ("\ngshhs_magics at %g km:\n# of points reduced by %.1f%% (out %d, in %d)\n# of polygons reduced by %.1f%% out (%d, in %d)\n", tolerance, redux, n_tot_out, n_tot_in, redux2, n_out, n_id);

      //exit (EXIT_SUCCESS); // C???
      return 0;
}

/*! \brief Stack-based Douglas Peucker line simplification routine
 returned value is the number of output points
 Kindly provided by  Dr. Gary J. Robinson,
 Environmental Systems Science Centre,
 University of Reading, Reading, UK
*/
int Douglas_Peucker_i (int x_source[], int y_source[], int n_source, double band, int index[])
/* x_source[]: Input coordinates in micro-degrees     */
/* y_source[]:                                  */
/* n_source:      Number of points              */
/* band:    tolerance in kilometers             */
/* index[]: output co-ordinates indices         */
{
      int   n_stack, n_dest, start, end, i, sig;
      int   *sig_start, *sig_end;   /* indices of start&end of working section */

      double dev_sqr, max_dev_sqr, band_sqr;
      double  x12, y12, d12, x13, y13, d13, x23, y23, d23;

        /* check for simple cases */

        if ( n_source < 3 ) return(0);    /* one or two points */

        /* more complex case. initialize stack */

      sig_start = (int *) get_memory (VNULL, n_source, sizeof (int), "Douglas_Peucker_i");
      sig_end   = (int *) get_memory (VNULL, n_source, sizeof (int), "Douglas_Peucker_i");
      
      band *= 360.0 / (2.0 * M_PI * 6371.007181);     /* Now in degrees */
      band_sqr = sqr(band);
            
      n_dest = 0;

        sig_start[0] = 0;
        sig_end[0] = n_source-1;

        n_stack = 1;

        /* while the stack is not empty  ... */

        while ( n_stack > 0 )
        {
                /* ... pop the top-most entries off the stacks */

                start = sig_start[n_stack-1];
                end = sig_end[n_stack-1];

                n_stack--;

                if ( end - start > 1 )  /* any intermediate points ? */
                {
                        /* ... yes, so find most deviant intermediate point to
                               either side of line joining start & end points */

                x12 = 1.0e-6 * (x_source[end] - x_source[start]);
                if (fabs (x12) > 180.0) x12 = 360.0 - fabs (x12);
                y12 = 1.0e-6 * (y_source[end] - y_source[start]);
            x12 *= cos (F * (y_source[end] + y_source[start]));
            d12 = sqr(x12) + sqr(y12);

                for ( i = start + 1, sig = start, max_dev_sqr = -1.0; i < end; i++ )
                        {
                                x13 = 1.0e-6 * (x_source[i] - x_source[start]);
                                y13 = 1.0e-6 * (y_source[i] - y_source[start]);
                        if (fabs (x13) > 180.0) x13 = 360.0 - fabs (x13);
                                x13 *= cos (F * (y_source[i] + y_source[start]));

                                x23 = 1.0e-6 * (x_source[i] - x_source[end]);
                                y23 = 1.0e-6 * (y_source[i] - y_source[end]);
                        if (fabs (x23) > 180.0) x23 = 360.0 - fabs (x23);
                                x23 *= cos (F * (y_source[i] + y_source[end]));
                                
                                d13 = sqr(x13) + sqr(y13);
                                d23 = sqr(x23) + sqr(y23);

                                if ( d13 >= ( d12 + d23 ) )
                                        dev_sqr = d23;
                                else if ( d23 >= ( d12 + d13 ) )
                                        dev_sqr = d13;
                                else
                                        dev_sqr =  sqr( x13 * y12 - y13 * x12 ) / d12;

                                if ( dev_sqr > max_dev_sqr  )
                                {
                                        sig = i;
                                        max_dev_sqr = dev_sqr;
                                }
                        }

                        if ( max_dev_sqr < band_sqr )   /* is there a sig. intermediate point ? */
                        {
                                /* ... no, so transfer current start point */

                                index[n_dest] = start;
                                n_dest++;
                        }
                        else
                        {
                                /* ... yes, so push two sub-sections on stack for further processing */

                                n_stack++;

                                sig_start[n_stack-1] = sig;
                                sig_end[n_stack-1] = end;

                                n_stack++;

                                sig_start[n_stack-1] = start;
                                sig_end[n_stack-1] = sig;
                        }
                }
                else
                {
                        /* ... no intermediate points, so transfer current start point */

                        index[n_dest] = start;
                        n_dest++;
                }
        }

        /* transfer last point */
        index[n_dest] = n_source-1;
        n_dest++;

      free ((void *)sig_start);
      free ((void *)sig_end);

        return (n_dest);
}

void *get_memory (void *prev_addr, int n, size_t size, char *progname)
{
        void *tmp;

      if (n == 0) return(VNULL); /* Take care of n = 0 */

      if (prev_addr) {
            if ((tmp = realloc ((void *) prev_addr, (size_t) (n * size))) == VNULL) {
                  fprintf (stderr, "gshhs Fatal Error: %s could not reallocate more memory, n = %d\n", progname, n);
                  exit (EXIT_FAILURE);
            }
      }
      else {
            if ((tmp = calloc ((size_t) n, (size_t) size)) == VNULL) {
                  fprintf (stderr, "gshhs Fatal Error: %s could not allocate memory, n = %d\n", progname, n);
                  exit (EXIT_FAILURE);
            }
      }
      return (tmp);
}

Generated by  Doxygen 1.6.0   Back to index