*** empty log message ***

Plugin/GSHHS.cpp

+// $Id: GSHHS.cpp,v 1.1 2008-03-18 09:27:59 remacle Exp $
+//
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// 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; either version 2 of the License, or
+// (at your option) any later version.
+//
+// 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.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+// 
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "GSHHS.h"
+#include "Field.h"
+#include "GModel.h"
+#include "GeoStringInterface.h"
+#include <sstream>
+#include <iostream>
+#include "Context.h"
+#include "GeoEarthImport.h"
+extern Context_T CTX;
+
+/*	$Id: GSHHS.cpp,v 1.1 2008-03-18 09:27:59 remacle Exp $
+ *
+ * PROGRAM:	gshhs.c
+ * AUTHOR:	Paul Wessel (pwessel@hawaii.edu)
+ * CREATED:	JAN. 28, 1996
+ * PURPOSE:	To extract ASCII data from binary shoreline data
+ *		as described in the 1996 Wessel & Smith JGR Data Analysis Note.
+ * VERSION:	1.1 (Byte flipping added)
+ *		1.2 18-MAY-1999:
+ *		   Explicit binary open for DOS systems
+ *		   POSIX.1 compliant
+ *		1.3 08-NOV-1999: Released under GNU GPL
+ *		1.4 05-SEPT-2000: Made a GMT supplement; FLIP no longer needed
+ *		1.5 14-SEPT-2004: Updated to deal with latest GSHHS database (1.3)
+ *
+ *	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 */
+
+/* For byte swapping on little-endian systems (GSHHS is bigendian) */
+#define swabi2(i2) (((i2) >> 8) + (((i2) & 255) << 8))
+#define swabi4(i4) (((i4) >> 24) + (((i4) >> 8) & 65280) + (((i4) & 65280) << 8) + (((i4) & 255) << 24))
+#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
+
+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 level;			/* 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake */
+	int west, east, south, north;	/* min/max extent in micro-degrees */
+	int area;			/* Area of polygon in 1/10 km^2 */
+	int version;			/* Version of GSHHS polygon (3 is latest and first with this item) */
+	short int greenwich;		/* Greenwich is 1 if Greenwich is crossed */
+	short int source;		/* 0 = CIA WDBII, 1 = WVS */
+};
+
+struct	POINT {	/* Each lon, lat pair is stored in micro-degrees in 4-byte integer format */
+	int	x;
+	int	y;
+};
+
+StringXNumber GSHHSOptions_Number[] = {
+  {GMSH_FULLRC, "iField", NULL, -1.}
+};
+StringXString GSHHSOptions_String[] = {
+	{GMSH_FULLRC, "InFileName",NULL,"gshhs_c.b"},
+	{GMSH_FULLRC, "OutFileName",NULL,"earth.geo"}
+};
+
+extern "C"
+{
+  GMSH_Plugin *GMSH_RegisterGSHHSPlugin()
+  {
+    return new GMSH_GSHHSPlugin();
+  }
+}
+
+void GMSH_GSHHSPlugin::getName(char *name) const
+{
+  strcpy(name, "GSHHS");
+}
+
+void GMSH_GSHHSPlugin::getInfos(char *author, char *copyright,
+				    char *help_text) const
+{
+  strcpy(author, "J. Lambrechts (jonathanlambrechts@gmail.org)");
+  strcpy(copyright, "GPL");
+  strcpy(help_text,
+	 "Plugin(GSHHS) import GSHHS data.\n");
+}
+
+int GMSH_GSHHSPlugin::getNbOptions() const
+{
+  return sizeof(GSHHSOptions_Number) / sizeof(StringXNumber);
+}
+int GMSH_GSHHSPlugin::getNbOptionsStr() const
+{
+  return sizeof(GSHHSOptions_String) / sizeof(StringXString);
+}
+
+StringXNumber *GMSH_GSHHSPlugin::getOption(int iopt)
+{
+  return &GSHHSOptions_Number[iopt];
+}
+StringXString *GMSH_GSHHSPlugin::getOptionStr(int iopt)
+{
+  return &GSHHSOptions_String[iopt];
+}
+
+void GMSH_GSHHSPlugin::catchErrorMessage(char *errorMessage) const
+{
+  strcpy(errorMessage, "GSHHS failed...");
+}
+
+PView *GMSH_GSHHSPlugin::execute(PView *v)
+{
+  int iField = (int)GSHHSOptions_Number[0].def;
+	char *filename=(char*)GSHHSOptions_String[0].def;
+	char *outfilename=(char*)GSHHSOptions_String[1].def;
+	Field *field=NULL;
+	GeoEarthImport geo_import(outfilename);
+	if(iField!=-1){
+		field = GModel::current()->fields.get(iField);
+		if(!field){
+			Msg(GERROR, "Field[%d] does not exist", iField);
+			return NULL;
+		}else{
+			geo_import.set_size_field(field);
+		}
+	}
+	FILE *fp;
+	if ((fp = fopen (filename, "rb")) == NULL ) {
+		Msg(GERROR, "gshhs: Could not find file %s.\n", filename);
+		return NULL;
+	}
+	double w, e, s, n;
+	char source;
+	int	k, max_east = 270000000, info, n_read, flip;
+	struct	POINT p;
+	struct GSHHS h;
+	while(1){
+		n_read = fread ((void *)&h, (size_t)sizeof (struct GSHHS), (size_t)1, fp);
+		if(n_read!=1)break;
+		flip = (! (h.level > 0 && h.level < 5));	/* Take as sign that byte-swabbing is needed */
+		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);
+		}
+		w = h.west  * 1.0e-6;	/* Convert from microdegrees to degrees */
+		e = h.east  * 1.0e-6;
+		s = h.south * 1.0e-6;
+		n = h.north * 1.0e-6;
+		source = (h.source == 1) ? 'W' : 'C';	/* Either WVS or CIA (WDBII) pedigree */
+		if ( h.level!=1) {	// Skip data (lake,island,pond)
+			fseek (fp, (long)(h.n * sizeof(struct POINT)), SEEK_CUR);
+			continue;
+		} 
+		for (k = 0; k < h.n; k++) {
+			if (fread ((void *)&p, (size_t)sizeof(struct POINT), (size_t)1, fp) != 1) {
+				Msg(GERROR, "gshhs:  Error reading gshhs file for polygon %d, point %d.\n",h.id, k);
+				return NULL;
+			}
+			if (flip) {
+				p.x = swabi4 ((unsigned int)p.x);
+				p.y = swabi4 ((unsigned int)p.y);
+			}
+			SPoint2 ll(
+				((h.greenwich && p.x > max_east) ? p.x * 1.0e-6 - 360.0 : p.x * 1.0e-6)*M_PI/180,
+				 (p.y * 1.0e-6)*M_PI/180);
+			geo_import.add_point_lon_lat(ll);
+		}
+		geo_import.end_loop();
+	} 
+	geo_import.end_surface();
+  return NULL;
+}

Plugin/GSHHS.h

+#ifndef _GSHHS_H_
+#define _GSHHS_H_
+
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// 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; either version 2 of the License, or
+// (at your option) any later version.
+//
+// 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.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+// 
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Plugin.h"
+
+extern "C"
+{
+  GMSH_Plugin *GMSH_RegisterGSHHSPlugin();
+}
+
+class GMSH_GSHHSPlugin : public GMSH_Post_Plugin
+{
+ public:
+  void getName(char *name) const;
+  void getInfos(char *author, char *copyright, char *help_text) const;
+  void catchErrorMessage(char *errorMessage) const;
+  int getNbOptions() const;
+  int getNbOptionsStr() const;
+  StringXNumber* getOption(int iopt);  
+  StringXString* getOptionStr(int iopt);  
+  PView *execute(PView *);
+};
+
+#endif

Plugin/GeoEarthImport.cpp

+#include "GeoEarthImport.h"
+GeoEarthImport::~GeoEarthImport(){
+	file.close();
+}
+GeoEarthImport::GeoEarthImport(const std::string _filename){
+	filename=_filename;
+	file.open(filename.c_str());
+	size_field=NULL;
+	loop_buff.precision(16);
+	std::ostringstream buff;
+	il=ip=ill=is=ifi=0;
+	buff<<"IP = newp;\n";
+	buff<<"IL = newl;\n";
+	buff<<"ILL = newll;\n";
+	buff<<"IS = news;\n";
+	buff<<"IFI = newf;\n";
+	buff<<"Point ( IP + "<<ip++<<" ) = {0, 0, 0 , 1};\n";
+	buff<<"Point ( IP + "<<ip++<<" ) = {0, 0, 6.371e6 , 1};\n";
+	buff<<"PolarSphere ( IS + "<<is++<<" ) = {IP , IP+1};\n";
+	file<<buff.str();	
+	new_surface();
+	new_loop();
+}
+void GeoEarthImport::set_size_field(Field *_size_field){
+	size_field=_size_field;
+}
+void GeoEarthImport::add_point(const SPoint3 &point){
+	SPoint3 midpoint = point;
+	SPoint2 stereo(-point.x()/(1+point.z()),-point.y()/(1+point.z()));
+	midpoint+=lastpoint;
+	midpoint/=2;
+	if( ip==first_point_in_loop || !size_field || point.distance(lastpoint)*6361e3 > (*size_field)(midpoint[0],midpoint[1],midpoint[2])){
+		loop_buff<<"Point ( IP + "<<ip++<<" ) = {"<<stereo.x()<<", "<<stereo.y()<<", "<<0<<", 1};\n";
+		lastpoint=point;
+	}
+}
+void GeoEarthImport::add_point_lon_lat(const SPoint2 &ll){
+	SPoint3 xyz( cos(ll.y())*cos(ll.x()), cos(ll.y())*sin(ll.x()), sin(ll.y()));
+	add_point(xyz);
+}
+void GeoEarthImport::new_loop(){
+	loop_buff.str("");
+	first_point_in_loop=ip;
+}
+void GeoEarthImport::end_loop(){
+	if(ip-first_point_in_loop>3){
+		loop_buff<<"BSpline ( IL + "<<il++<<" ) = { IP + "<<first_point_in_loop<<" : IP + "<<ip-1<<", IP + "<<first_point_in_loop<<" };\n";
+		loop_buff<<"Line Loop ( ILL + "<<ill++<<" ) = { IL + "<<il-1<<" };";
+		file<<loop_buff.str();	
+		if(!empty_surface)
+			surface_buff<<", ";
+		surface_buff<< "ILL + "<<ill-1;
+		empty_surface=false;
+	}else{
+		ip=first_point_in_loop;
+	}
+	new_loop();
+}
+void GeoEarthImport::new_surface(){
+	surface_buff.str("");
+	surface_buff<<"Plane Surface( IS + "<<is++<<" ) = { ";
+	first_point_in_surface=ip;
+	bool empty_surface=true;
+}
+void GeoEarthImport::end_surface(){
+	if(!empty_surface){
+		surface_buff<<"};\n";
+		surface_buff<<"Field [ IFI + "<<ifi<<"]  = Attractor;\n";
+		surface_buff<<"Field [ IFI + "<<ifi<<"].NodesList  = { IP + "<<first_point_in_surface <<" : IP + "<<ip-1<<" };";
+		ifi++;
+		file<<surface_buff.str();	
+	}
+	new_surface();
+}
+
+/*static void projectLatLon(Point3D &psphere,Point &pplan){
+	pplan.x=asin(psphere.z);
+	pplan.y=atan2(psphere.y,psphere.x);
+	pplan.zone_id=psphere.zone_id;
+}
+static void projectInvStereo(Point &pplan,Point3D &psphere){
+	double rp2=pplan.x*pplan.x+pplan.y*pplan.y;
+	psphere.z= (4-rp2)/(4+rp2);
+	psphere.x= (1+psphere.z)/2*pplan.x;
+	psphere.y= (1+psphere.z)/2*pplan.y;
+	psphere.zone_id=pplan.zone_id;
+}
+static void projectStereo(Point3D &psphere,Point &pplan){
+	pplan.x=2*psphere.x/(psphere.z+1);
+	pplan.y=2*psphere.y/(psphere.z+1);
+	pplan.zone_id=psphere.zone_id;
+}*/

Plugin/GeoEarthImport.h

+#ifndef _GEO_EARTH_IMPORT_H_
+#define _GEO_EARTH_IMPORT_H_
+#include "SPoint3.h"
+#include "Field.h"
+#include <sstream>
+#include <fstream>
+#include "GeoStringInterface.h"
+class GeoEarthImport{
+	Field *size_field;
+	std::ostringstream loop_buff,surface_buff;;
+	std::string filename;
+	std::ofstream file;
+	int il,ip,is,ill,ifi;
+	int first_point_in_loop,first_point_in_surface;
+	bool empty_surface;
+	SPoint3 lastpoint;
+	void new_surface();
+	void new_loop();
+
+	public :
+	GeoEarthImport(const std::string _filename);
+	~GeoEarthImport();
+	void set_size_field(Field *_size_field);
+	void add_point(const SPoint3 &point);
+	void add_point_lon_lat(const SPoint2 &ll);
+	void end_loop();
+	void end_surface();
+};
+#endif