[GRASS-SVN] r35643 - grass-addons/vector/v.in.gshhs
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jan 26 06:42:45 EST 2009
Author: mmetz
Date: 2009-01-26 06:42:45 -0500 (Mon, 26 Jan 2009)
New Revision: 35643
Added:
grass-addons/vector/v.in.gshhs/Makefile
grass-addons/vector/v.in.gshhs/README
grass-addons/vector/v.in.gshhs/description.html
grass-addons/vector/v.in.gshhs/gshhs.h
grass-addons/vector/v.in.gshhs/gshhs_datum.txt
grass-addons/vector/v.in.gshhs/main.c
Log:
upload GSHHS import module
Added: grass-addons/vector/v.in.gshhs/Makefile
===================================================================
--- grass-addons/vector/v.in.gshhs/Makefile (rev 0)
+++ grass-addons/vector/v.in.gshhs/Makefile 2009-01-26 11:42:45 UTC (rev 35643)
@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.in.gshhs
+
+LIBES = $(VECTLIB) $(GPROJLIB) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(VECTDEP) $(GPROJDEP) $(DBMIDEP) $(GISDEP)
+
+EXTRA_INC = $(VECT_INC) $(PROJINC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Added: grass-addons/vector/v.in.gshhs/README
===================================================================
--- grass-addons/vector/v.in.gshhs/README (rev 0)
+++ grass-addons/vector/v.in.gshhs/README 2009-01-26 11:42:45 UTC (rev 35643)
@@ -0,0 +1,47 @@
+
+######################
+v.in.gshhs
+
+Program to extract and import binned shoreline data
+available from:
+
+http://www.soest.hawaii.edu/wessel/gshhs/gshhs.html
+
+The Global Self Consistant Hierarchical High Resolution Shorelines
+(GSHHS) are available at 5 resolutions.
+
+- crude
+- low
+- intermediate
+- high
+- full
+
+Bob Covill
+
+######################
+
+5/2005:
+
+ A new version of GSHHS was published in 9/2004 with a slightly
+ modified file format. The v.in.gshhs was updated accordingly.
+
+Markus Neteler
+
+######################
+
+1/2009:
+
+ v.in.gshhs has been updated to work with grass-6.4
+ this module has been tested with GSHHS versions 6, 10, and 11
+
+ further modifications:
+
+ the bounding box of shorelines is now used to check if a line should be
+ imported
+
+ a small attribute table allowing fast queries is written for layer 1
+ indicating the shoreline type
+ category values for layer 2 are set to GSHHS line id
+
+Markus Metz
+
Added: grass-addons/vector/v.in.gshhs/description.html
===================================================================
--- grass-addons/vector/v.in.gshhs/description.html (rev 0)
+++ grass-addons/vector/v.in.gshhs/description.html 2009-01-26 11:42:45 UTC (rev 35643)
@@ -0,0 +1,71 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.in.gshhs</em> imports GSHHS shorelines as available on
+<a href="ftp://ftp.soest.hawaii.edu/pwessel/gshhs/">ftp://ftp.soest.hawaii.edu/pwessel/gshhs/</a>
+<p>
+GSHHS shorelines are in latlon but can be imported to any location with
+any projection apart from unreferenced XY projection. Reprojection is
+done on the fly during import.
+<p>
+Shorelines are imported as lines, not boundaries.
+<p>
+The categories in layer 1 refer to the type of shoreline: shore, land,
+lake, island in lake, pond in island in lake. The field <tt>type</tt> is
+set to these descriptions. All lines of the same type have the same
+category. The categories in layer 2 refer to the GSHHS ID of each
+imported line. Each line has a unique category value (the GSHHS ID) in
+layer 2, that may be useful for further processing. An attribute table
+for layer 2 is not created.
+<p>
+The import region is restricted to either the current computational
+region or to the extend specified with the options n, s, e, w. Either of
+the two possibilities must be used. All lines that overlap with the
+specified region are imported, similar to <em>v.in.ogr</em>.
+
+<h2>NOTES</h2>
+The GSHHS shorelines are in files named <tt>gshhs_[f|h|i|l|c].b</tt> where
+the letter in brackets indicates the resolution. Only these files can be
+imported, import of world borders and world rivers also included in the zip
+archive is currently not supported.
+<p>
+The GSHHS shorelines are available at 5 resolutions:
+<ul>
+<li>f: full resolution
+<li>h: high resolution
+<li>i: intermediate resolution
+<li>l: low resolution
+<li>c: coarse resolution
+</ul>
+<p>
+The generated table for layer 1 allows a fast query of the shoreline type.
+If needed a table for layer 2 can be added with <em>v.db.addtable</em>.
+The new table can be populated with category values from layer 2 (unique
+line ID) with <em>v.to.db</em>. Shoreline type can be uploaded from
+layer 1 to layer 2 with <em>v.to.db</em>.
+<p>
+Long lines can be split with <em>v.split</em>.
+<p>
+Lines can be converted to boundaries with <em>v.type</em>.
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a HREF="v.db.addtable.html">v.db.addtable</a>,
+<a HREF="v.split.html">v.split</a>,
+<a HREF="v.to.db.html">v.to.db</a>,
+<a HREF="v.type.html">v.type</a>,
+</em>
+
+<h2>REFERENCES</h2>
+
+The processing and assembly of the GSHHS data is described in
+<p>
+Wessel, P., and W. H. F. Smith, 1996. A Global Self-consistent, Hierarchical,
+High-resolution Shoreline Database. J. Geophys. Res., 101(B4), 8741-8743.
+
+<h2>AUTHORS</h2>
+The gshhstograss tool was written by Simon Cox and Paul Wessel.<br>
+The original version of v.in.gshhs was written by Bob Covill based on gshhstograss.<br>
+Modifications and updates by Markus Neteler and Markus Metz.
+<p>
+<i>Last changed: $Date: 2009-01-26 18:43:35 +0200 (Mon, 26 Jan 2009) $</i>
Added: grass-addons/vector/v.in.gshhs/gshhs.h
===================================================================
--- grass-addons/vector/v.in.gshhs/gshhs.h (rev 0)
+++ grass-addons/vector/v.in.gshhs/gshhs.h 2009-01-26 11:42:45 UTC (rev 35643)
@@ -0,0 +1,75 @@
+/* $Id: gshhs.h,v 1.17 2008/01/23 03:22:49 guru Exp $
+ *
+ * Include file defining structures used in gshhs.c
+ *
+ * Paul Wessel, SOEST
+ *
+ * Copyright (c) 1996-2008 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
+ *
+ * 14-SEP-2004. PW: Version 1.3. Header is now n * 8 bytes (n = 5)
+ * For use with version 1.3 of GSHHS
+ * 2-MAY-2006. PW: Version 1.4. Header is now 32 bytes (all int 4)
+ * For use with version 1.4 of GSHHS
+ * 31-MAR-2007. PW: Version 1.5. no format change
+ * For use with version 1.5 of GSHHS
+ * 28-AUG-2007. PW: Version 1.6. no format change
+ * For use with version 1.6 of GSHHS which now has WDBII
+ * borders and rivers.
+ */
+
+#ifndef _GSHHS
+#define _GSHHS
+#define _POSIX_SOURCE 1 /* GSHHS code is POSIX compliant */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+#ifndef SEEK_CUR /* For really ancient systems */
+#define SEEK_CUR 1
+#endif
+
+#define GSHHS_DATA_VERSION 6 /* For v1.5 data set */
+#define GSHHS_PROG_VERSION "1.9"
+
+#define GSHHS_SCL 1.0e-6 /* COnvert micro-degrees to degrees */
+
+/* For byte swapping on little-endian systems (GSHHS is defined to be bigendian) */
+
+#define swabi4(i4) (((i4) >> 24) + (((i4) >> 8) & 65280) + (((i4) & 65280) << 8) + (((i4) & 255) << 24))
+
+struct GSHHS { /* Global Self-consistent Hierarchical High-resolution Shorelines */
+ int id; /* Unique polygon id number, starting at 0 */
+ int n; /* Number of points in this polygon */
+ int flag; /* = level + version << 8 + greenwich << 16 + source << 24 */
+ /* flag contains 4 items, one in each byte, as follows:
+ * low byte: level = flag & 255: Values: 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake
+ * 2nd byte: version = (flag >> 8) & 255: Values: Should be 4 for GSHHS version 1.4
+ * 3rd byte: greenwich = (flag >> 16) & 255: Values: Greenwich is 1 if Greenwich is crossed
+ * 4th byte: source = (flag >> 24) & 255: Values: 0 = CIA WDBII, 1 = WVS
+ */
+ int west, east, south, north; /* min/max extent in micro-degrees */
+ int area; /* Area of polygon in 1/10 km^2 */
+};
+
+struct GSHHS_POINT { /* Each lon, lat pair is stored in micro-degrees in 4-byte integer format */
+ int x;
+ int y;
+};
+#endif /* _GSHHS */
Property changes on: grass-addons/vector/v.in.gshhs/gshhs.h
___________________________________________________________________
Name: svn:executable
+
Added: grass-addons/vector/v.in.gshhs/gshhs_datum.txt
===================================================================
--- grass-addons/vector/v.in.gshhs/gshhs_datum.txt (rev 0)
+++ grass-addons/vector/v.in.gshhs/gshhs_datum.txt 2009-01-26 11:42:45 UTC (rev 35643)
@@ -0,0 +1,97 @@
+Leming, T.D., May, L.N., Jones, P. 1999. A Geographic Information System for
+Near Reat-Time use of Remote Sensing in Fisheries Management in the Gulf of
+Mexico. 85 p. February 11, 1999.
+http://coastwatch.noaa.gov/pubs/digdocs/nrtgis_final_report2.pdf
+
+"A continental land mask was created from the global, self-consistent,
+hierarchical, high-resolution shoreline database (GSHHS) compiled by Wessel
+and Smith (1996). The GSHHS is public domain digital coastline file available
+from the School of Ocean and Earth Science Technology at the University of
+Hawaii12 and the NOAA National Geophysical Data Center 13 (NGDC) and contains
+coastlines and closed waterbodies extracted from two widely used public domain
+databases, the World Data Bank II coastline (WDB) and the World Vector
+Shoreline (WVS). Lohrenz (1988) indicated that the working scale of the WVS
+data is about about 1:250,000 and is referenced to the World Geodetic System
+of 1972 (WGS-72) horizontal datum and the mean high water (MHW) vertical
+datum."
+
+Lohrenz, M.C. 1988. Design and implementation of the digital world vector
+shoreline data format. Report 194, Naval Ocean Research and Development
+Activity, Stennis Space Center, Miss., 34 pp.
+
+#######################################################################
+http://wegener.mechanik.tu-darmstadt.de/GMT-Help/Archiv/1960.html
+Coastlines, etc
+
+ * This message: [ Message body ] [ More options ]
+ * Related messages: [ Next message ] [ Previous message ]
+
+From: Paul Wessel <pal.wessel_at_>
+Date: 1999-04-14 08:11:02
+
+Hi GMT users-
+
+Each year I get a few questions like these:
+
+1. "I am making a plot of a small area, using pscoast -Df, and the island that
+ I used to visit as a kid does not show up (or lake, or river, etc)"
+
+2. "When I compared the coastline from pscoast with GPS measurements, the
+ coastline seems shifted by a few hundred meters"
+
+3. "Where this river is wide it is missing entirely, with the result that
+ some of the islands in the river show up as lakes"
+
+To simplify the answer, I will assume you have read Appendix K in the GMT
+Reference & Cookbook section on how the GSHHS coastlines were put together.
+
+Answers:
+
+ There are many possibilities:
+ a) The feature was missing in the original WVS (land outlines) or WDBII
+ (lakes, rivers, boundaries) data set.
+ b) The feature was eaten by our preprocessing steps.
+ c) WVS and WDBII were sometimes incompatible.
+ d) Different datums were used
+
+ One must remember that the WVS coastline is much higher accuracy than the
+ WDBII lakes and rivers. It has happened that a feature that was considered
+ a coastal lake in WVS (i.e., it did not break the shoreline and thus was
+ not included in WVS) was a bay in WDBII (and thus not included as a lake).
+ The result is that the lake never appeared in GSHHS.
+ Likewise, rivers do not always end at the coastline; it may cross it or not
+ quite make it to the ocean.
+
+How can we improve the data? Ideally, we would prefer that some agency
+announces a sucessor to WVS and WDBII that is better, and we can start from
+that point. We do not want to be the curators of these data. In the mean time:
+
+1. If you can show us that a missing feature is indeed in the original
+ WVS or WDBII data then _maybe_ I can find time to check the processing
+code.
+ I need much motivation to dig into all that stuff again.
+
+2. If you have a high-resolution (i.e., comparable to other features in -Df)
+ digital version of the missing feature, please email it to us with an
+ explanatory note. We may then incorporate it in the next version.
+
+3. The rivers seem to be messed up, either in the source or by us. If you
+ have digitized missing pieces please share them with us.
+
+4. The shift that some have documented may be related to datum. If anybody
+ knows what reference ellipsoid was used for WVS and WDBII we would like
+ to hear it. We made no correction to the lon/lat when making the data
+ compilation.
+
+5. We still miss some recent political borders. It may be premature to
+ add lines in the Balkans, but if you have high-quality political borders
+ not in GMT, please share with us.
+
+Paul Wessel, Professor
+Dept. of Geology & Geophysics
+School of Ocean & Earth Science & Technology
+University of Hawaii at Manoa
+1680 East-West Road,
+Honolulu, HI 96822
+(808) 956-4778/5154 (voice/fax)
+
Added: grass-addons/vector/v.in.gshhs/main.c
===================================================================
--- grass-addons/vector/v.in.gshhs/main.c (rev 0)
+++ grass-addons/vector/v.in.gshhs/main.c 2009-01-26 11:42:45 UTC (rev 35643)
@@ -0,0 +1,432 @@
+/**********************************************************************
+ *
+ * MODULE: v.in.gshhs (based on gshhstograss.c)
+ *
+ * AUTHORS: Simon Cox (simon at ned.dem.csiro.au) & (gshhstograss.c)
+ * Paul Wessel (wessel at soest.hawaii.edu) & (gshhstograss.c)
+ * Bob Covill <bcovill at tekmap.ns.ca> (v.in.gshhs)
+ * Markus Metz <markus.metz.giswork googlemail.com> (v.in.gshhs for grass 6.4)
+ *
+ * PURPOSE: To extract GRASS binary vector data from binary gshhs
+ * shoreline data as described in
+ * Wessel, P., and W.H.F. Smith. A Global Self-consistent,
+ * Hierarchical, High-resolution Shoreline Database,
+ * J. Geophys. Res., 101(B4), 8741-8743, 1996.
+ * http://www.soest.hawaii.edu/wessel/hshhs/gshhs.html
+ *
+ * LICENSE: This program is free software under the
+ * GNU General Public License (>=v2).
+ * Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ **********************************************************************/
+
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <time.h>
+
+#include <grass/gis.h>
+#include <grass/dbmi.h>
+#include <grass/Vect.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+/* updating to a newer GSHHS version can be as easy as replacing gshhs.h
+ * if not too many changes were introduced */
+#include "gshhs.h"
+
+int main(int argc, char **argv)
+{
+ double w, e, s, n, area, lon, lat;
+ double minx = -180., maxx = 180., miny = -90., maxy = 90.;
+ char *dataname, *outname;
+ char buf[2000], source;
+ static char *slevel[] =
+ { "shore", "land", "lake", "island in lake", "pond in island in lake" };
+ int shore_levels = 5;
+ FILE *fp;
+ int k, i, max_east = 270000000, n_read, flip, level, version, greenwich, src;
+ int cnt = 1;
+ struct GSHHS_POINT p; /* renamed to avoid conflict */
+ struct GSHHS h;
+ struct pj_info info_in;
+ struct pj_info info_out;
+ struct Key_Value *out_proj_keys, *out_unit_keys;
+ int type, zone;
+ struct Cell_head region;
+ struct Map_info VectMap;
+ struct line_pnts *Points;
+ struct line_cats *Cats;
+ double lat_nw, lon_nw, lat_ne, lon_ne;
+ double lat_sw, lon_sw, lat_se, lon_se;
+ struct
+ {
+ struct Option *input, *output, *n, *s, *e, *w;
+ } parm;
+ struct
+ {
+ struct Flag *g, *topo, *a;
+ } flag;
+ struct GModule *module;
+
+ /* attribute table */
+ struct field_info *Fi;
+ dbDriver *driver;
+ dbString sql, strval;
+ char *cat_col_name = "cat";
+
+
+ G_gisinit(argv[0]);
+
+ /* Set description */
+ module = G_define_module();
+ module->description = ""
+ "Imports Global Self-consistent Hierarchical High-resolution Shoreline (GSHHS) vector data";
+
+ parm.input = G_define_option();
+ parm.input->key = "input";
+ parm.input->type = TYPE_STRING;
+ parm.input->required = YES;
+ parm.input->gisprompt = "old_file,file,gshhs";
+ parm.input->description = "Name of GSHHS shoreline file: gshhs_[f|h|i|l|c].b";
+
+ parm.output = G_define_standard_option(G_OPT_V_OUTPUT);
+
+ parm.n = G_define_option();
+ parm.n->key = "n";
+ parm.n->type = TYPE_DOUBLE;
+ parm.n->required = NO;
+ parm.n->description = "North limit in units of location projection";
+
+ parm.s = G_define_option();
+ parm.s->key = "s";
+ parm.s->type = TYPE_DOUBLE;
+ parm.s->required = NO;
+ parm.s->description = "South limit in units of location projection";
+
+ parm.e = G_define_option();
+ parm.e->key = "e";
+ parm.e->type = TYPE_DOUBLE;
+ parm.e->required = NO;
+ parm.e->description = "East limit in units of location projection";
+
+ parm.w = G_define_option();
+ parm.w->key = "w";
+ parm.w->type = TYPE_DOUBLE;
+ parm.w->required = NO;
+ parm.w->description = "West limit in units of location projection";
+
+ flag.g = G_define_flag();
+ flag.g->key = 'g';
+ flag.g->description = "Get coordinates from current GRASS region";
+
+ flag.topo = G_define_flag();
+ flag.topo->key = 't';
+ flag.topo->description =
+ "Do not build topology for output vector";
+
+/* Importing as boundary causes problems with subregions, incorrect boundaries */
+/* flag.a = G_define_flag();
+ flag.a->key = 'a';
+ flag.a->description =
+ "Import shoreline as type line (default = boundary)";
+*/
+
+ if (G_parser(argc, argv))
+ exit(1);
+
+/* get parameters */
+ dataname = parm.input->answer;
+ outname = parm.output->answer;
+
+ G_get_window(®ion);
+ zone = region.zone;
+
+/* Out Info */
+ out_proj_keys = G_get_projinfo();
+ out_unit_keys = G_get_projunits();
+ if (pj_get_kv(&info_out, out_proj_keys, out_unit_keys) < 0) {
+ exit(0);
+ }
+ G_free_key_value(out_proj_keys);
+ G_free_key_value(out_unit_keys);
+
+/* In Info */
+/* N.B. GSHHS data is not referenced to any ellipsoid or datum. This
+ * limits its precision to several hundred metres. Hence it does not
+ * matter which ellipsoid we use in the input projection keys as the
+ * precision of the data is less than the maximum error that could be
+ * caused by projecting using the 'wrong' ellipsoid.
+ * Anyone who understands this better than me feel free to change it.
+ * PK March 2003 */
+
+/* GSHHS spatial reference
+ * GSHHS is a combination of WVS and WDBII
+ * The WVS dataset uses WGS84 according to
+ * http://shoreline.noaa.gov/data/datasheets/wvs.html
+ * The WDBII dataset may have used wgs72 but is inaccurate anyway
+ * Safest would be to always use wgs84
+ * MM Jan 2009 */
+
+/* set input projection to lat/long w/ same ellipsoid as output */
+/* TODO: use wgs84 datum and ellipsoid */
+ info_in.zone = 0;
+ info_in.meters = 1.;
+ sprintf(info_in.proj, "ll");
+ if ((info_in.pj = pj_latlong_from_proj(info_out.pj)) == NULL)
+ G_fatal_error("Unable to set up lat/long projection parameters");
+
+ if (flag.g->answer) {
+ /* get coordinates from current region */
+ minx = region.west;
+ maxx = region.east;
+ miny = region.south;
+ maxy = region.north;
+ }
+ else {
+ /* get coordinates from command line */
+ if (!parm.w->answer || !parm.e->answer || !parm.s->answer
+ || !parm.n->answer) {
+ G_fatal_error("Missing region parameters: you must supply n, s, e, w values or use '-g' flag");
+ }
+
+ sscanf(parm.w->answer, "%lf", &minx);
+ sscanf(parm.e->answer, "%lf", &maxx);
+ sscanf(parm.s->answer, "%lf", &miny);
+ sscanf(parm.n->answer, "%lf", &maxy);
+
+ if (maxx < minx)
+ G_fatal_error("East must be east of West");
+
+ if (maxy < miny)
+ G_fatal_error("North must be north of South");
+ }
+
+
+ if (G_projection() != PROJECTION_LL) {
+ /* convert to latlon and print bounds */
+
+ /* NW */
+ lon_nw = minx;
+ lat_nw = maxy;
+ if (pj_do_proj(&lon_nw, &lat_nw, &info_out, &info_in) < 0) {
+ G_fatal_error("Error in coordinate transformation for north west corner");
+ }
+ /* NE */
+ lon_ne = maxx;
+ lat_ne = maxy;
+ if (pj_do_proj(&lon_ne, &lat_ne, &info_out, &info_in) < 0) {
+ G_fatal_error("Error in coordinate transformation for north east corner");
+ }
+ /* SE */
+ lon_se = maxx;
+ lat_se = miny;
+ if (pj_do_proj(&lon_se, &lat_se, &info_out, &info_in) < 0) {
+ G_fatal_error("Error in coordinate transformation for south eest corner");
+ }
+ /* SW */
+ lon_sw = minx;
+ lat_sw = miny;
+ if (pj_do_proj(&lon_sw, &lat_sw, &info_out, &info_in) < 0) {
+ G_fatal_error("Error in coordinate transformation for south west corner");
+ }
+
+ /* get north max */
+ maxy = lat_nw > lat_ne ? lat_nw : lat_ne;
+ /* get south min*/
+ miny = lat_sw < lat_se ? lat_sw : lat_se;
+ /* get east max*/
+ maxx = lon_ne > lon_se ? lon_ne : lon_se;
+ /* get west min*/
+ minx = lon_nw < lon_sw ? lon_nw : lon_sw;
+ }
+
+ /* check coordinate accuracy */
+ if (maxx < minx) {
+ G_fatal_error(_("maxx %f < minx %f."), maxx, minx);
+ }
+
+ if (maxy < miny) {
+ G_fatal_error(_("maxy %f < miny %f."), maxy, miny);
+ }
+
+ G_message(_("Using lat/lon bounds N=%f S=%f E=%f W=%f\n"), maxy, miny, maxx,
+ minx);
+
+ /* open GSHHS shoreline for reading */
+ if ((fp = fopen(dataname, "rb")) == NULL) {
+ G_fatal_error(_("Could not find file %s."), dataname);
+ }
+
+ /* Open new vector */
+ if (0 > Vect_open_new(&VectMap, outname, 0)) {
+ G_fatal_error(_("Cannot open new vector %s"), outname);
+ }
+
+ /* set vector line type to GV_LINE, boundaries can cause problems */
+ type = GV_LINE;
+
+ /* Initialize vector line struct */
+ Points = Vect_new_line_struct();
+
+ /* Initialize vector category struct */
+ Cats = Vect_new_cats_struct();
+
+ /* read header from GSHHS database */
+ n_read = fread ((void *)&h, (size_t)sizeof (struct GSHHS), (size_t)1, fp);
+ version = (h.flag >> 8) & 255;
+ flip = (version < 6); /* Take as sign that byte-swapping is needed */
+ if (flip) {
+ version = swabi4 ((unsigned int)h.flag);
+ version = (version >> 8) & 255;
+ }
+
+ /* read lines from GSHHS database */
+ while (n_read == 1) {
+ if (flip) {
+ h.id = swabi4 ((unsigned int)h.id);
+ h.n = swabi4 ((unsigned int)h.n);
+ 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.flag = swabi4 ((unsigned int)h.flag);
+ }
+ level = h.flag & 255;
+ /* version = (h.flag >> 8) & 255; */
+ greenwich = (h.flag >> 16) & 255;
+ src = (h.flag >> 24) & 255;
+ w = h.west * GSHHS_SCL; /* Convert from microdegrees to degrees */
+ e = h.east * GSHHS_SCL;
+ s = h.south * GSHHS_SCL;
+ n = h.north * GSHHS_SCL;
+ /* skip WDBII data because they are too inacurate? */
+ source = (src == 1) ? 'W' : 'C'; /* Either WVS or CIA (WDBII) pedigree */
+ area = 0.1 * h.area; /* Now im km^2 */
+
+ if (h.west > max_east) w -= 360.0;
+ if (h.east > max_east) e -= 360.0;
+
+ /* got line header now read in points if line box overlaps with region */
+ if( ( w <= maxx && e >= minx ) && ( s <= maxy && n >= miny ) ){
+ Vect_reset_line(Points);
+ Vect_reset_cats(Cats);
+
+ for (k = 0; k < h.n; k++) {
+ if (fread
+ ((void *)&p, (size_t)sizeof(struct GSHHS_POINT), (size_t)1,
+ fp) != 1)
+ G_fatal_error("Error reading data.");
+
+ if (flip) {
+ p.x = swabi4((unsigned int) p.x);
+ p.y = swabi4((unsigned int) p.y);
+ }
+ lon = p.x * GSHHS_SCL;
+ if (p.x > max_east) lon -= 360.0;
+ lat = p.y * GSHHS_SCL;
+
+ if (G_projection() != PROJECTION_LL) {
+ if (pj_do_proj(&lon, &lat, &info_in, &info_out) < 0) {
+ G_fatal_error("Error in coordinate transformation");
+ }
+ }
+ Vect_append_point(Points, lon, lat, 0.);
+ } /* done with line */
+ if (k) {
+ /* simple table and cats for layer 1: not unique, just per slevel */
+ Vect_cat_set(Cats, 1, level);
+ /* cats for layer 2 are unique ID, same like GSHHS ID */
+ Vect_cat_set(Cats, 2, h.id);
+ Vect_write_line(&VectMap, type, Points, Cats);
+ cnt++;
+ }
+ }
+ else {
+ G_debug(1, "skipped line box west: %f, east: %f, north: %f, south: %f", w, e, n, s);
+ fseek (fp, (long)(h.n * sizeof(struct GSHHS_POINT)), SEEK_CUR);
+ }
+ max_east = 180000000; /* Only Eurasiafrica needs 270 */
+
+ n_read = fread((void *)&h, (size_t)sizeof (struct GSHHS), (size_t)1, fp);
+ }
+ /* done with gshhs file */
+ fclose(fp);
+
+ Vect_destroy_line_struct(Points);
+ Vect_destroy_cats_struct(Cats);
+
+ /* write out vect header info */
+ Vect_set_scale(&VectMap, 100000);
+ sprintf(buf, "GSHHS shorelines version %d imported by v.in.gshhs", version);
+ Vect_set_comment(&VectMap, buf);
+
+ /* create table and link new vector to table, code taken from v.in.ogr */
+
+ db_init_string(&sql);
+ db_init_string(&strval);
+
+ Fi = Vect_default_field_info(&VectMap, 1, NULL, GV_1TABLE);
+
+ Vect_map_add_dblink(&VectMap, 1, NULL, Fi->table,
+ cat_col_name, Fi->database, Fi->driver);
+
+ sprintf(buf, "create table %s (%s integer, type varchar(25))", Fi->table,
+ cat_col_name);
+ db_set_string(&sql, buf);
+
+ driver = db_start_driver_open_database(Fi->driver,
+ Vect_subst_var(Fi->database, &VectMap));
+
+ if (driver == NULL) {
+ G_fatal_error(_("Cannot open database %s by driver %s"),
+ Vect_subst_var(Fi->database, &VectMap), Fi->driver);
+ }
+
+ G_debug(1, "table: %s", Fi->table);
+ G_debug(1, "driver: %s", Fi->driver);
+ G_debug(1, "database: %s", Fi->database);
+
+ if (db_execute_immediate(driver, &sql) != DB_OK) {
+ db_close_database(driver);
+ db_shutdown_driver(driver);
+ G_fatal_error(_("Cannot create table: %s"),
+ db_get_string(&sql));
+ }
+
+ if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
+ G_warning(_("Cannot create index"));
+
+ if (db_grant_on_table(driver, Fi->table, DB_PRIV_SELECT,
+ DB_GROUP | DB_PUBLIC) != DB_OK)
+ G_fatal_error(_("Cannot grant privileges on table %s"), Fi->table);
+
+ db_begin_transaction(driver);
+
+ for (i = 0; i < shore_levels; i++) {
+ sprintf(buf, "insert into %s values ( %d, \'%s\')", Fi->table, i + 1, slevel[i]);
+ db_set_string(&sql, buf);
+
+ if (db_execute_immediate(driver, &sql) != DB_OK) {
+ db_close_database(driver);
+ db_shutdown_driver(driver);
+ G_fatal_error(_("Cannot insert new row: %s"),
+ db_get_string(&sql));
+ }
+ }
+ db_commit_transaction(driver);
+ db_close_database_shutdown_driver(driver);
+
+ if (flag.topo->answer) {
+ G_message(_("Run \"v.build map=%s option=build\" to build topology"), outname);
+ }
+ else {
+ Vect_build(&VectMap);
+ }
+ Vect_close(&VectMap);
+
+ G_message(_("GSHHS shoreline version %d import complete"), version);
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/vector/v.in.gshhs/main.c
___________________________________________________________________
Name: svn:executable
+
More information about the grass-commit
mailing list