Skip to content
Snippets Groups Projects
Commit f6fcc770 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

-add an option in GSHHS Plugin to remove narrow straits (almost complete rewrite of GSHHS plugin) 

-complete the doc of the GSHHS Plugin
-compile GSHHS plugin even if matheval is not present
parent 2c102dff
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Plugin/GSHHS.cpp
+ 907
− 439
View file @ f6fcc770
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include "GSHHS.h"
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "Field.h"
#include "GModel.h"
#include <sstream>
#include <fstream>
/*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;
}*/
// ************** Coordinate Systems *************
class coordinate_system{
public:
virtual void to_cartesian(const SPoint3 coord, SPoint3 &cartesian)=0;
virtual void from_cartesian(const SPoint3 cartesian, SPoint3 &coord)=0;
};
// ************** Longitude Latitude ***************
class coordinate_lonlat:public coordinate_system{
public:
void to_cartesian(const SPoint3 ll, SPoint3 &xyz){
double clat = cos(ll.y());
xyz.setPosition( clat * cos(ll.x()), clat * sin(ll.x()), sin(ll.y()));
}
void from_cartesian(const SPoint3 xyz, SPoint3 &ll){
double r=sqrt(xyz.x()*xyz.x()+xyz.y()*xyz.y()+xyz.z()*xyz.z());
ll.setPosition( atan2(xyz.y(),xyz.x()), asin(xyz.z()/r), 0);
}
};
// ************** Longitude Latitude (degrees) ***************
class coordinate_lonlat_degrees:public coordinate_system{
coordinate_lonlat cll;
SPoint3 llradian;
public:
void to_cartesian(const SPoint3 ll, SPoint3 &xyz){
llradian.setPosition(ll.x()*M_PI/180,ll.y()*M_PI/180,0);
cll.to_cartesian(llradian,xyz);
}
void from_cartesian(const SPoint3 xyz, SPoint3 &ll){
cll.from_cartesian(xyz,llradian);
ll.setPosition(llradian.x()*180/M_PI,llradian.y()*180/M_PI,0);
}
};
#include "SPoint2.h"
#include "GSHHS.h"
#include "GModel.h"
// ************** UTM **************
class coordinate_utm :public coordinate_system{
int zone;
coordinate_lonlat ll_conv;
SPoint3 ll;
double a, b, n, n2, n3, n4, n5, e, e2, e1, e12, e13, e14, J1, J2, J3, J4,
Ap, Bp, Cp, Dp, Ep, e4, e6, ep, ep2, ep4, k0, mu_fact;
public:
static int get_zone_from_longitude(double lon){
return (int)ceil((lon / M_PI + 1) * 30);
}
double meridionalarc(double lon, double lat){
return Ap * lat + Bp * sin(2 * lat) + Cp * sin(4 * lat) +
Dp * sin(6 * lat) + Ep;
}
void from_cartesian(const SPoint3 xyz,SPoint3 &utm) {
ll_conv.from_cartesian(xyz,ll);
double S = meridionalarc(ll.x(),ll.y());
double slat = sin(ll.y());
double clat = cos(ll.y());
double slat2 = slat * slat;
double clat2 = clat * clat;
double clat3 = clat2 * clat;
double clat4 = clat3 * clat;
double tlat2 = slat2 / clat2;
double nu = a / sqrt(1 - e * e * slat2);
double p = ll.x() - ((zone - 0.5) / 30 - 1) * M_PI;
double p2 = p * p;
double p3 = p * p2;
double p4 = p2 * p2;
utm.setPosition(
k0 * nu * clat * p + (k0 * nu * clat3 / 6)
* (1 - tlat2 + ep2 * clat2) * p3 + 5e5,
S * k0 + k0 * nu * slat * clat / 2 * p2
+ k0 * nu * slat * clat3 / 24
* (5 - tlat2 + 9 * ep2 * clat2 + 4 * ep4 * clat4) * p4,
0);
}
void to_cartesian(const SPoint3 utm, SPoint3 &xyz){
double mu = utm.y() * mu_fact;
double fp =
mu + J1 * sin(2 * mu) + J2 * sin(4 * mu) + J3 * sin(6 * mu) +
J4 * sin(8 * mu);
double cfp = cos(fp);
double cfp2 = cfp * cfp;
double sfp = sin(fp);
double sfp2 = sfp * sfp;
double c1 = ep2 * cfp2;
double c12 = c1 * c1;
double t1 = sfp2 / cfp2;
double t12 = t1 * t1;
double r1 = a * (1 - e2) / pow((1 - e2 * sfp2), 1.5);
double n1 = a / sqrt(1 - e2 * sfp2);
double d = (utm.x() - 5e5) / (n1 * k0);
double d2 = d * d;
double d3 = d2 * d;
double d4 = d2 * d2;
double d5 = d4 * d;
double d6 = d4 * d2;
ll.setPosition(
((zone - 0.5) / 30 - 1) * M_PI + (d - (1 + 2 * t1 + c1) * d3 / 6 +
(5 - 2 * c1 + 28 * t1 - 3 * c12 +
8 * ep2 +
24 * t12) * d5 / 120) / cfp,
fp - n1 * sfp / cfp / r1
* (d2 / 2 - (5 + 3 * t1 + 10 * c1 - 4 * c12 - 9 * ep2) * d4 / 24
+ (61 + 90 * t1 + 298 * c1 + 45 * t12 - 3 * c12 - 252 * ep2) * d6 / 720),
0);
ll_conv.to_cartesian(ll,xyz);
}
coordinate_utm(int _zone, double _a = 6378137, double _b = 6356752.3142) {
/* see http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.HTM */
a = _a; /* Equatorial Radius*/
b = _b; /* Rayon Polar Radius*/
zone=_zone;
n = (a - b) / (a + b);
n2 = n * n;
n3 = n * n * n;
n4 = n * n * n * n;
n5 = n * n * n * n * n;
e = sqrt(1 - b * b / a / a);
e2 = e * e;
e1 = (1 - sqrt(1 - e2)) / (1 + sqrt(1 - e2));
e12 = e1 * e1;
e13 = e1 * e1 * e1;
e14 = e1 * e1 * e1 * e1;
J1 = (3 * e1 / 2 - 27 * e13 / 32);
J2 = (21 * e12 / 16 - 55 * e14 / 32);
J3 = 151 * e13 / 96;
J4 = 1097 * e14 / 512;
Ap = a * (1 - n + (5. / 4.) * (n2 - n3) + (81. / 64.) * (n4 - n5));
Bp = -3 * a * n / 2 * (1 - n + (7. / 8.) * (n2 - n3) +
(55. / 64.) * (n4 - n5));
Cp = 14 * a * n2 / 16 * (1 - n + (3. / 4) * (n2 - n3));
Dp = -35 * a * n3 / 48 * (1 - n + 11. / 16. * (n2 - n3));
Ep = +315 * a * n4 / 51 * (1 - n);
e4 = e2 * e2;
e6 = e2 * e2 * e2;
ep = e * a / b;
ep2 = ep * ep;
ep4 = ep2 * ep2;
k0 = 0.9996;
mu_fact = 1 / (k0 * a * (1 - e2 / 4 - 3 * e4 / 64 - 5 * e6 / 256));
}
};
class GeoEarthImport
class GMSH_GSHHSPlugin:public GMSH_Post_Plugin
{
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,first_point_in_attractor;
bool empty_surface,closed_loop;
SPoint3 lastpoint;
coordinate_system *c_syst;
void new_attractor(){
first_point_in_attractor = ip;
}
void new_surface() {
surface_buff.str("");
surface_buff << "Plane Surface( IS + " << is++ << " ) = { ";
first_point_in_surface = ip;
empty_surface = true;
}
void new_loop(){
loop_buff.str("");
first_point_in_loop = ip;
closed_loop=true;
}
public:
void set_coordinate_system(coordinate_system *_c_syst){
c_syst=_c_syst;
}
GeoEarthImport(const std::string _filename, bool write_polar_sphere){
c_syst=NULL;
filename = _filename;
file=new std::ofstream(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";
if(write_polar_sphere){
buff << "Point ( IP + " << ip++ << " ) = {0, 0, 0 };\n";
buff << "Point ( IP + " << ip++ <<" ) = {0, 0, 6.371e6 };\n";
buff << "PolarSphere ( IS + " << is++ << " ) = {IP , IP+1};\n";
}
*file << buff.str();
new_surface();
new_attractor();
new_loop();
}
~GeoEarthImport(){
//file << "Euclidian Coordinates;";
file->close();
}
void set_size_field(Field * _size_field){
size_field = _size_field;
}
void add_point(SPoint3 point){
if(c_syst)
c_syst->to_cartesian(point,point);
double lc=0;
if(size_field){
double d=(*size_field)(point[0]* 6361e3,point[1]* 6361e3,point[2]* 6361e3);
if(d<0){
if(first_point_in_loop!=ip){
closed_loop=false;
end_loop();
closed_loop=false;
// ************** Inputs (readers) *************
class reader{
public:
virtual int next_loop(bool &closed)=0;
virtual bool next_point(SPoint3 &point)=0;
};
class reader_loops2:public reader{
FILE *fp;
int npoints_in_loop;
int ipoint;
std::string filename;
public:
reader_loops2(std::string _filename){
filename=_filename;
fp=fopen(filename.c_str(),"r");
}
~reader_loops2(){
fclose(fp);
}
int next_loop(bool & closed){
ipoint=0;
npoints_in_loop=-1;
int i_closed;
if(fscanf(fp, "%d %d", &npoints_in_loop,&i_closed) != 2)
return 0;
closed = i_closed;
return npoints_in_loop;
}
bool next_point(SPoint3 &point){
if(ipoint>=npoints_in_loop)
return false;
point[2]=0;
ipoint++;
if(fscanf(fp, "%le %le", &point[0], &point[1]) != 2)
Msg::Error("gshhs: Error reading loops2 file.\n");
return true;
}
};
class reader_gshhs:public reader{
/* $Id: GSHHS.cpp,v 1.19 2009-01-12 11:59:39 remacle 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.
*/
static const int GSHHS_DATA_VERSION=6; // For v1.5 data set
//define GSHHS_PROG_VERSION "1.9"
static const double GSHHS_SCL=1.0e-6; // COnvert micro-degrees to degrees
inline unsigned int swabi4(unsigned int i4){ // For byte swapping on little-endian systems (GSHHS is defined to be bigendian)
return (((i4) >> 24) + (((i4) >> 8) & 65280) + (((i4) & 65280) << 8) + (((i4) & 255) << 24));
}
class POINT { /* Each lon, lat pair is stored in micro-degrees in 4-byte integer format */
public:
int x;
int y;
};
class GSHHS { /* Global Self-consistent Hierarchical High-resolution Shorelines */
public:
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 */
};
GSHHS h;
POINT p;
FILE *fp;
int max_east,flip,ip,greenwich;
bool first_loop;
public:
reader_gshhs(std::string filename){
fp=fopen(filename.c_str(),"r");
max_east = 270000000;
first_loop=true;
}
~reader_gshhs(){
fclose(fp);
}
int next_loop(bool &closed){
closed=true;
int level=0;
ip=0;
while(level!=1){
int n_read = fread ((void *)&h, (size_t)sizeof (struct GSHHS), (size_t)1, fp);
if(n_read!=1 || feof(fp))
return 0;
level = h.flag & 255;
int version = (h.flag >> 8) & 255;
flip = (version != GSHHS_DATA_VERSION); /* 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.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);
}
if(level!=1)
fseek(fp, (size_t)(h.n * sizeof(struct POINT)), SEEK_CUR);
if(first_loop)
first_loop=false;
else
max_east = 180000000; /* Only Eurasiafrica needs 270 */
}
return;
greenwich = (h.flag >> 16) & 255;
int src = (h.flag >> 24) & 255;
double w = h.west * GSHHS_SCL; /* Convert from microdegrees to degrees */
double e = h.east * GSHHS_SCL;
double s = h.south * GSHHS_SCL;
double n = h.north * GSHHS_SCL;
char source = (src == 1) ? 'W' : 'C'; /* Either WVS or CIA (WDBII) pedigree */
int line = (h.area) ? 0 : 1; /* Either Polygon (0) or Line (1) (if no area) */
double area = 0.1 * h.area; /* Now im km^2 */
return h.n;
}
bool next_point(SPoint3 &point){
if(ip>=h.n)
return false;
if (fread ((void *)&p, (size_t)sizeof(struct POINT), (size_t)1, fp) != 1) {
printf ("gshhs: Error reading gshhs file.\n");
exit(1);
}
if (flip) {
p.x = swabi4 ((unsigned int)p.x);
p.y = swabi4 ((unsigned int)p.y);
}
double lon = p.x * GSHHS_SCL;
if (greenwich && p.x > max_east)
lon -= 360.0;
double lat = p.y * GSHHS_SCL;
point[0]=lon*M_PI/180;
point[1]=lat*M_PI/180;
point[2]=0;
ip++;
return true;
}
};
// ************** Coordinate Systems *************
class coordinate_system{
public:
virtual void to_cartesian(const SPoint3 coord, SPoint3 &cartesian)=0;
virtual void from_cartesian(const SPoint3 cartesian, SPoint3 &coord)=0;
};
// ************** Longitude Latitude ***************
class coordinate_lonlat:public coordinate_system{
double radius;
public:
coordinate_lonlat(double r){
radius=r;
}
void to_cartesian(const SPoint3 ll, SPoint3 &xyz){
double clat = cos(ll.y());
xyz.setPosition( clat * cos(ll.x())*radius, clat * sin(ll.x())*radius, sin(ll.y())*radius);
}
void from_cartesian(const SPoint3 xyz, SPoint3 &ll){
double r=sqrt(xyz.x()*xyz.x()+xyz.y()*xyz.y()+xyz.z()*xyz.z());
ll.setPosition( atan2(xyz.y(),xyz.x()), asin(xyz.z()/r), r);
}
};
// ************** Longitude Latitude (degrees) ***************
class coordinate_lonlat_degrees:public coordinate_system{
coordinate_lonlat cll;
SPoint3 llradian;
public:
coordinate_lonlat_degrees(double r):cll(r){};
void to_cartesian(const SPoint3 ll, SPoint3 &xyz){
llradian.setPosition(ll.x()*M_PI/180,ll.y()*M_PI/180,0);
cll.to_cartesian(llradian,xyz);
}
lc=std::min(d, (*size_field) (lastpoint[0],lastpoint[1],lastpoint[2]));
void from_cartesian(const SPoint3 xyz, SPoint3 &ll){
cll.from_cartesian(xyz,llradian);
ll.setPosition(llradian.x()*180/M_PI,llradian.y()*180/M_PI,0);
}
};
// ************** UTM **************
class coordinate_utm :public coordinate_system{
int zone;
coordinate_lonlat ll_conv;
SPoint3 ll;
double a, b, n, n2, n3, n4, n5, e, e2, e1, e12, e13, e14, J1, J2, J3, J4,
Ap, Bp, Cp, Dp, Ep, e4, e6, ep, ep2, ep4, k0, mu_fact;
public:
static int get_zone_from_longitude(double lon){
return (int)ceil((lon / M_PI + 1) * 30);
}
double meridionalarc(double lon, double lat){
return Ap * lat + Bp * sin(2 * lat) + Cp * sin(4 * lat) +
Dp * sin(6 * lat) + Ep;
}
void from_cartesian(const SPoint3 xyz,SPoint3 &utm) {
ll_conv.from_cartesian(xyz,ll);
double S = meridionalarc(ll.x(),ll.y());
double slat = sin(ll.y());
double clat = cos(ll.y());
double slat2 = slat * slat;
double clat2 = clat * clat;
double clat3 = clat2 * clat;
double clat4 = clat3 * clat;
double tlat2 = slat2 / clat2;
double nu = a / sqrt(1 - e * e * slat2);
double p = ll.x() - ((zone - 0.5) / 30 - 1) * M_PI;
double p2 = p * p;
double p3 = p * p2;
double p4 = p2 * p2;
utm.setPosition(
k0 * nu * clat * p + (k0 * nu * clat3 / 6)
* (1 - tlat2 + ep2 * clat2) * p3 + 5e5,
S * k0 + k0 * nu * slat * clat / 2 * p2
+ k0 * nu * slat * clat3 / 24
* (5 - tlat2 + 9 * ep2 * clat2 + 4 * ep4 * clat4) * p4,
0);
}
void to_cartesian(const SPoint3 utm, SPoint3 &xyz){
double mu = utm.y() * mu_fact;
double fp =
mu + J1 * sin(2 * mu) + J2 * sin(4 * mu) + J3 * sin(6 * mu) +
J4 * sin(8 * mu);
double cfp = cos(fp);
double cfp2 = cfp * cfp;
double sfp = sin(fp);
double sfp2 = sfp * sfp;
double c1 = ep2 * cfp2;
double c12 = c1 * c1;
double t1 = sfp2 / cfp2;
double t12 = t1 * t1;
double r1 = a * (1 - e2) / pow((1 - e2 * sfp2), 1.5);
double n1 = a / sqrt(1 - e2 * sfp2);
double d = (utm.x() - 5e5) / (n1 * k0);
double d2 = d * d;
double d3 = d2 * d;
double d4 = d2 * d2;
double d5 = d4 * d;
double d6 = d4 * d2;
ll.setPosition(
((zone - 0.5) / 30 - 1) * M_PI + (d - (1 + 2 * t1 + c1) * d3 / 6 +
(5 - 2 * c1 + 28 * t1 - 3 * c12 +
8 * ep2 +
24 * t12) * d5 / 120) / cfp,
fp - n1 * sfp / cfp / r1
* (d2 / 2 - (5 + 3 * t1 + 10 * c1 - 4 * c12 - 9 * ep2) * d4 / 24
+ (61 + 90 * t1 + 298 * c1 + 45 * t12 - 3 * c12 - 252 * ep2) * d6 / 720),
0);
ll_conv.to_cartesian(ll,xyz);
}
coordinate_utm(int _zone,double ll_radius, double _a = 6378137, double _b = 6356752.3142) :ll_conv(ll_radius){
/* see http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.HTM */
a = _a; /* Equatorial Radius*/
b = _b; /* Rayon Polar Radius*/
zone=_zone;
n = (a - b) / (a + b);
n2 = n * n;
n3 = n * n * n;
n4 = n * n * n * n;
n5 = n * n * n * n * n;
e = sqrt(1 - b * b / a / a);
e2 = e * e;
e1 = (1 - sqrt(1 - e2)) / (1 + sqrt(1 - e2));
e12 = e1 * e1;
e13 = e1 * e1 * e1;
e14 = e1 * e1 * e1 * e1;
J1 = (3 * e1 / 2 - 27 * e13 / 32);
J2 = (21 * e12 / 16 - 55 * e14 / 32);
J3 = 151 * e13 / 96;
J4 = 1097 * e14 / 512;
Ap = a * (1 - n + (5. / 4.) * (n2 - n3) + (81. / 64.) * (n4 - n5));
Bp = -3 * a * n / 2 * (1 - n + (7. / 8.) * (n2 - n3) +
(55. / 64.) * (n4 - n5));
Cp = 14 * a * n2 / 16 * (1 - n + (3. / 4) * (n2 - n3));
Dp = -35 * a * n3 / 48 * (1 - n + 11. / 16. * (n2 - n3));
Ep = +315 * a * n4 / 51 * (1 - n);
e4 = e2 * e2;
e6 = e2 * e2 * e2;
ep = e * a / b;
ep2 = ep * ep;
ep4 = ep2 * ep2;
k0 = 0.9996;
mu_fact = 1 / (k0 * a * (1 - e2 / 4 - 3 * e4 / 64 - 5 * e6 / 256));
}
};
/********** classes and functions to ensure minimal distance and angle between points **********/
class box;
class loop;
class point{
public:
SPoint3 v;
std::list<point>::iterator it_loop;
std::list<point*>::iterator it_box;
box *b;
double min_dist;
loop *l;
point(double _x,double _y,double _z,Field *f,double straits_factor){
v[0]=_x; v[1]=_y; v[2]=_z;
if(f)
min_dist=(*f)(v[0],v[1],v[2])*straits_factor;
}
point(double _x,double _y,double _z,double _min_dist){
v[0]=_x; v[1]=_y; v[2]=_z;
min_dist=_min_dist;
}
double dist(point p){
return sqrt(
(v[0]-p.v[0])*(v[0]-p.v[0])
+(v[1]-p.v[1])*(v[1]-p.v[1])
+(v[2]-p.v[2])*(v[2]-p.v[2]));
}
void to_stereo(double &xp,double &yp,bool inverse_stereo=false){
if(inverse_stereo){
yp=v[1]/(1-v[2]);
xp=v[0]/(1-v[2]);
}else{
xp=-v[1]/(1+v[2]);
yp=-v[0]/(1+v[2]);
}
}
void to_latlon(double &lat,double &lon){
lat=asin(v[2]);
lon=atan2(v[1],v[0]);
if(lon<-15*M_PI/16)lon+=2*M_PI;
}
};
class box{
public:
box *sub_boxes[2];
point min,max;
int max_n;
double min_size,size;
bool splitted;
std::list<point*> list;
int cut_dim,ndim;
int n;
box(point _min,point _max,int _max_n=100,double _min_size=0,int _cut_dim=0):min(_min),max(_max),max_n(_max_n){
n=0;
min_size=_min_size;
size=0;
ndim=3;
for(int i=0;i<ndim;i++)
size=std::max(size,fabs(max.v[i]-min.v[i]));
splitted=false;
cut_dim=_cut_dim;
}
void add( point &p){
if(!splitted){
list.push_back(&p);
p.it_box=list.end();
p.it_box--;
p.b=this;
n++;
if(n>max_n && size>=min_size){
point mid0(min);
point mid1(max);
mid1.v[cut_dim]=mid0.v[cut_dim]=(min.v[cut_dim]+max.v[cut_dim])/2;
int newdim=(cut_dim+1)%ndim;
sub_boxes[0]=new box(min,mid1,max_n,min_size,newdim);
sub_boxes[1]=new box(mid0,max,max_n,min_size,newdim);
splitted=true;
for(std::list<point*>::iterator it=list.begin();it!=list.end();it++)
add(**it);
list.clear();
}
}else{
int ix=(p.v[cut_dim]-min.v[cut_dim])>(max.v[cut_dim]-min.v[cut_dim])/2;
sub_boxes[ix]->add(p);
}
}
void remove(point *p){
if(p->b!=this)
p->b->remove(p);
else{
list.erase(p->it_box);
n--;
}
}
void find_closest(point p,int nc,double *d,point **cp){
int i=(cut_dim+ndim-1)%ndim;
if(p.v[i]<min.v[i]-d[nc-1] || p.v[i]>max.v[i]+d[nc-1])
return;
if(splitted){
for(int i=0;i<2;i++)
sub_boxes[i]->find_closest(p,nc,d,cp);
}else for(std::list<point*>::iterator it=list.begin();it!=list.end();it++){
double dd=p.dist(**it);
if(dd<d[nc-1]){
int j;
for(j=nc-2;j>=0;j--){
if(dd>d[j])break;
d[j+1]=d[j];
cp[j+1]=cp[j];
}
d[j+1]=dd;
cp[j+1]=*it;
}
}
}
};
static double get_angle(double x0,double y0,double x1,double y1,double x2,double y2){
double dx0=x1-x0;
double dx1=x2-x1;
double dy0=y1-y0;
double dy1=y2-y1;
return atan2(dx0*dy1-dx1*dy0,dx1*dx0+dy1*dy0);
}
lc/=6361e3;
if(ip == first_point_in_loop || point.distance(lastpoint) >lc){
SPoint2 stereo(-point.x() / (1 + point.z()), -point.y() / (1 + point.z()));
loop_buff << "Point ( IP + " << ip++ << " ) = {" << stereo.
x() << ", " << stereo.y() << ", " << 0 << " };\n";
lastpoint = point;
static double get_angle(point &p0,point &p1,point &p2){
double dx0=p1.v[0]-p0.v[0];
double dy0=p1.v[1]-p0.v[1];
double dz0=p1.v[2]-p0.v[2];
double dx1=p2.v[0]-p1.v[0];
double dy1=p2.v[1]-p1.v[1];
double dz1=p2.v[2]-p1.v[2];
double scalar=dx0*dx1+dy0*dy1+dz0*dz1;
double vectx=dy0*dz1-dy1*dz0;
double vecty=dz0*dx1-dz1*dx0;
double vectz=dx0*dy1-dx1*dy0;
int sign=vectx*p1.v[0]+vecty*p1.v[1]+vectz*p1.v[2]<0?-1:1;
double vect=sqrt(vectx*vectx+vecty*vecty+vectz*vectz);
return atan2(-vect*sign,scalar);
}
}
void end_loop(bool closed=true) {
closed&=closed_loop;
if(ip - first_point_in_loop > 3) {
loop_buff<<"LoopStart"<<il<<" = IP + "<< first_point_in_loop<<";\n";
loop_buff<<"LoopEnd"<<il<<" = IP + "<< ip - 1<<";\n";
loop_buff << "BSpline ( IL + " << il++ << " ) = { IP + " <<
first_point_in_loop << " : IP + " << ip - 1 ;
if(closed) loop_buff<< ", IP + " << first_point_in_loop;
loop_buff<< " };\n";
if(closed){
loop_buff << "Line Loop ( ILL + " << ill++ << " ) = { IL + " << il -
1 << " };";
}
*file << loop_buff.str();
if(closed){
if(!empty_surface)
surface_buff << ", ";
surface_buff << "ILL + " << ill - 1;
empty_surface = false;
}
// angle (01,12)
class loop:public std::list<point>{
public:
bool closed;
loop(){
closed=true;
}
inline iterator next(iterator i){
++i;
if(i==end()){
i=begin();
}
return i;
}
inline iterator prev(iterator i){
if(i==begin())
i=end();
--i;
return i;
}
iterator remove_range(loop::iterator i0,loop::iterator i1){ //remove [id0, id1] and replace it by a single point at (id0+id1)/2;
for(int i=0;i<3;i++)
i1->v[i]=(i0->v[i]+i1->v[i])/2;
while(i0!=i1){
if(i0->b)
i0->b->remove(&*i0);
i0=erase(i0);
if(i0==end())
i0=begin();
}
return i1;
}
int orientation(iterator i0, iterator i1,bool reverse_stereo){
if(next(i0)==i1)
return 0;
double alpha=0;
double x[3],y[3];
i1->to_stereo(x[0],y[0],reverse_stereo);
i0->to_stereo(x[1],y[1],reverse_stereo);
iterator p=i0;
do{
p=next(p);
p->to_stereo(x[2],y[2],reverse_stereo);
alpha+=get_angle(x[0],y[0],x[1],y[1],x[2],y[2]);
x[0]=x[1];y[0]=y[1];
x[1]=x[2];y[1]=y[2];
}while(p!=i1);
i0->to_stereo(x[2],y[2],reverse_stereo);
alpha+=get_angle(x[0],y[0],x[1],y[1],x[2],y[2]);
return rint(alpha/(M_PI*2));
}
int length(iterator i0,iterator i1){
int l=1;
while(i0!=i1){
i0=next(i0);
l++;
}
return l;
}
void insert_range(iterator i,iterator j0,iterator j1){
iterator n=j0;
do{
j0=n;
iterator is=insert(i,*j0);
is->it_loop=is;
is->l=this;
*(j0->it_box)=&*is;
n=j0->l->next(j0);
}while(j0!=j1);
}
};
class loops:public std::list<loop>{
public:
void print_gnuplot(std::ostream &stream){
double lat,lon;
for(iterator il=begin();il!=end();il++){
stream<<"\n";
for(loop::iterator ip=il->begin();ip!=il->end();ip++){
ip->to_latlon(lat,lon);
stream<<lon*180/M_PI<<" "<<lat*180/M_PI<<"\n";
}
il->front().to_latlon(lat,lon);
stream<<lon*180/M_PI<<" "<<lat*180/M_PI<<"\n";
}
}
};
double stereo_cross_product(point p00,point p01, point p10,point p11){
double x[4],y[4];
p00.to_stereo(x[0],y[0]);
p01.to_stereo(x[1],y[1]);
p10.to_stereo(x[2],y[2]);
p11.to_stereo(x[3],y[3]);
return (x[1]-x[0])*(y[3]-y[2])-(y[1]-y[0])*(x[3]-x[2]);
}
else {
ip = first_point_in_loop;
bool is_intersected(point p00,point p01, point p10,point p11){
if(stereo_cross_product(p00,p01,p00,p10)*stereo_cross_product(p00,p01,p00,p11)>0)
return false;
if(stereo_cross_product(p10,p11,p10,p00)*stereo_cross_product(p10,p11,p10,p01)>0)
return false;
return true;
}
new_loop();
}
void end_surface(){
if(!empty_surface) {
surface_buff << "};\n";
surface_buff.str("");
*file << surface_buff.str()<<"\n";
void loop_fill_box(loop *l,box &b){
int i=0;
for(loop::iterator ip=l->begin();ip!=l->end();ip++){
ip->l=l;
ip->it_loop=ip;
b.add(*ip);
}
}
bool loop_check_intersections(loop *l,box &b){
// Check for intersections
bool result=false;
for(loop::iterator i00=l->begin();i00!=l->end();++i00){
loop::iterator i01=l->next(i00);
if(i01==l->begin() && !l->closed)
break;
double length0=i00->dist(*i01);
#define NP 5
double d[NP*2]={length0*1.0001};
for(int i=0;i<NP*2;i++)
d[i]=DBL_MAX;
point *cp[NP*2];
b.find_closest(*i00,NP,d,cp);
b.find_closest(*i01,NP,d+NP,cp+NP);
for(int i=0;i<NP*2;i++){
if(d[i]>length0)
break;
if(i00!=cp[i]->it_loop && i01!=cp[i]->it_loop){
loop::iterator i10=cp[i]->it_loop;
loop::iterator i11=l->next(i10);
if((i11!=l->begin() || l->closed) &&(i11 != i00 && i11 != i01))
if(is_intersected(*i00,*i01,*i10,*i11)){
if(l->length(i11,i00)<l->length(i01,i10))
i00=l->remove_range(i11,i00);
else
i00=l->remove_range(i01,i10);
result=true;
if(i00!=l->begin())
i00--;
break;
}
i11=i10;
i10=l->prev(i11);
if((i11!=l->begin() || l->closed) &&(i10 != i00 && i10 != i01))
if(is_intersected(*i00,*i01,*i10,*i11)){
if(l->length(i11,i00)<l->length(i01,i10))
i00=l->remove_range(i11,i00);
else
i00=l->remove_range(i01,i10);
result=true;
if(i00!=l->begin())
i00--;
break;
}
}
}
}
return result;
}
new_surface();
}
void end_attractor(){
*file << "Field [ IFI + " << ifi << "] = Attractor;\n";
*file << "Field [ IFI + " << ifi++ << "].NodesList = { IP + " <<
first_point_in_attractor << " : IP + " << ip - 1 << " };";
new_attractor();
}
};
/* $Id: GSHHS.cpp,v 1.18 2008-12-29 09:43:39 geuzaine 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.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
bool loop_check_close_points_self(loop *l,box &b){
bool result=false;
bool reverse_stereo=(l->orientation(l->begin(),--l->end(),false)<0);
for(loop::iterator i=l->begin();i!=l->end();){
double d[2]={i->min_dist*1.001,i->min_dist*1.001};
point *cp[2];
b.find_closest(*i,2,d,cp);
if(d[1]<i->min_dist){
loop::iterator id1=cp[1]->it_loop;
int not_a_loop=-1;
if(!l->closed){
loop::iterator ii;
for(ii=i;ii!=l->end() && ii!=id1;ii++);
not_a_loop= ii==l->end()?0:1;
}
if(not_a_loop!=0 && (l->orientation(i,id1,reverse_stereo)==-1 || l->length(i,id1)<3)){
i=l->remove_range(i,id1);
result=true;
}
if(not_a_loop!=1 &&(l->orientation(id1,i,reverse_stereo)==-1 || l->length(id1,i)<3)){
i=l->remove_range(id1,i);
result=true;
}else
i++;
}else
i++;
}
return result;
}
#define GSHHS_DATA_VERSION 6 /* For v1.5 data set */
#define GSHHS_PROG_VERSION "1.9"
bool loop_check_small_angles(loop *l){
bool removed=false;
for(loop::iterator i=l->begin();i!=l->end();++i){
loop::iterator i1=l->next(i);
loop::iterator i2=l->next(i1);
if((!l->closed) && i2==l->begin())
break;
double alpha=get_angle(*i,*i1,*i2);
if(alpha>3*M_PI/4){
i1->b->remove(&*i1);
i=l->erase(i1);
removed=true;
}
}
return removed;
}
#define GSHHS_SCL 1.0e-6 /* COnvert micro-degrees to degrees */
bool loop_check_close_points(loop *l,box &b){
#define NPD 10
for(loop::iterator i=l->begin();i!=l->end();){
double d[NPD];
point *cp[NPD];
for(int j=0;j<NPD;j++){
d[j]=i->min_dist*1.001;
cp[j]=NULL;
}
b.find_closest(*i,NPD,d,cp);
bool merged=false;
for(int j=0;j<NPD;j++){
if(!cp[j])break;
if(i->l!=cp[j]->l){
double lat,lon;
i->to_latlon(lat,lon);
loop::iterator f0=i;
loop::iterator f1=cp[j]->it_loop;
double newx[3],dx[3],crossx[3];
for(int k=0;k<3;k++){
newx[k]=(f0->v[k]+f1->v[k])/2;
dx[k]=(f1->v[k]-f0->v[k])/2;
}
for(int k=0;k<3;k++){
int k1=(k+1)%3;
int k2=(k+2)%3;
crossx[k]=newx[k1]*dx[k2]-newx[k2]*dx[k1];
}
double norm=sqrt(crossx[0]*crossx[0]+crossx[1]*crossx[1]+crossx[2]*crossx[2]);
for(int k=0;k<3;k++){
crossx[k]*=i->min_dist*0.1/norm;
f0->v[k]=newx[k]-crossx[k];
f1->v[k]=newx[k]+crossx[k];
}
/* For byte swapping on little-endian systems (GSHHS is defined to be bigendian) */
l->insert_range(++f0,f1->l->next(f1),f1);
cp[j]->l->clear();
merged=true;
break;
}
}
if(!merged)
i++;
}
}
class GeoEarthImport
{
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,first_point_in_attractor;
bool empty_surface;
void new_attractor(){
first_point_in_attractor = ip;
}
void new_surface() {
surface_buff.str("");
surface_buff << "Plane Surface( IS + " << is++ << " ) = { ";
first_point_in_surface = ip;
empty_surface = true;
}
void new_loop(){
loop_buff.str("");
first_point_in_loop = ip;
}
coordinate_lonlat cll;
public:
GeoEarthImport(const std::string _filename, bool write_polar_sphere,double radius):
cll(radius)
{
filename = _filename;
file=new std::ofstream(filename.c_str());
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";
if(write_polar_sphere){
buff << "Point ( IP + " << ip++ << " ) = {0, 0, 0 };\n";
buff << "Point ( IP + " << ip++ <<" ) = {0, 0,"<<radius<<"};\n";
buff << "PolarSphere ( IS + " << is++ << " ) = {IP , IP+1};\n";
}
*file << buff.str();
new_surface();
new_attractor();
new_loop();
}
~GeoEarthImport(){
//file << "Euclidian Coordinates;";
file->close();
}
void add_point(SPoint3 point){
SPoint3 pll;
cll.from_cartesian(point,pll);
SPoint2 stereo(-pll.x() / (1 + pll.z()), -pll.y() / (1 + pll.z()));
loop_buff << "Point ( IP + " << ip++ << " ) = {" << stereo.
x() << ", " << stereo.y() << ", " << 0 << " };\n";
}
void end_loop(bool closed) {
if(ip - first_point_in_loop > 3) {
loop_buff<<"LoopStart"<<il<<" = IP + "<< first_point_in_loop<<";\n";
loop_buff<<"LoopEnd"<<il<<" = IP + "<< ip - 1<<";\n";
loop_buff << "Spline ( IL + " << il++ << " ) = { IP + " <<
first_point_in_loop << " : IP + " << ip - 1 ;
if(closed) loop_buff<< ", IP + " << first_point_in_loop;
loop_buff<< " };\n";
if(closed){
loop_buff << "Line Loop ( ILL + " << ill++ << " ) = { IL + " << il -
1 << " };";
}
*file << loop_buff.str();
if(closed){
if(!empty_surface)
surface_buff << ", ";
surface_buff << "ILL + " << ill - 1;
empty_surface = false;
}
}
else {
ip = first_point_in_loop;
}
new_loop();
}
void end_surface(){
if(!empty_surface) {
surface_buff << "};\n";
surface_buff.str("");
*file << surface_buff.str()<<"\n";
}
new_surface();
}
void end_attractor(){
*file << "Field [ IFI + " << ifi << "] = Attractor;\n";
*file << "Field [ IFI + " << ifi++ << "].NodesList = { IP + " <<
first_point_in_attractor << " : IP + " << ip - 1 << " };";
new_attractor();
}
};
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 *);
};
#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 POINT { /* Each lon, lat pair is stored in micro-degrees in 4-byte integer format */
int x;
int y;
};
void import_gshhs(FILE * fp, GeoEarthImport & geo_import)
{
double w, e, s, n;
char source;
int k, max_east = 270000000, n_read, flip;
struct POINT p;
struct GSHHS h;
n_read = fread ((void *)&h, (size_t)sizeof (struct GSHHS), (size_t)1, fp);
int version = (h.flag >> 8) & 255;
flip = (version != GSHHS_DATA_VERSION); /* Take as sign that byte-swabbing is needed */
while (n_read == 1 && ! feof(fp)) {
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);
}
int level = h.flag & 255;
version = (h.flag >> 8) & 255;
int greenwich = (h.flag >> 16) & 255;
int 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;
source = (src == 1) ? 'W' : 'C'; /* Either WVS or CIA (WDBII) pedigree */
int line = (h.area) ? 0 : 1; /* Either Polygon (0) or Line (1) (if no area) */
double area = 0.1 * h.area; /* Now im km^2 */
if(level != 1) { // Skip data (lake,island,pond)
fseek(fp, (size_t)(h.n * sizeof(struct POINT)), SEEK_CUR);
}else{
for (k = 0; k < h.n; k++) {
if (fread ((void *)&p, (size_t)sizeof(struct POINT), (size_t)1, fp) != 1) {
Msg::Error ("gshhs: Error reading gshhs file.\n");
return;
}
if (flip) {
p.x = swabi4 ((unsigned int)p.x);
p.y = swabi4 ((unsigned int)p.y);
}
double lon = p.x * GSHHS_SCL;
if (greenwich && p.x > max_east) lon -= 360.0;
double lat = p.y * GSHHS_SCL;
geo_import.add_point(SPoint3(M_PI/180*lon,M_PI/180*lat,0));
}
geo_import.end_loop();
}
max_east = 180000000; /* Only Eurasiafrica needs 270 */
n_read = fread((void *)&h, (size_t)sizeof (struct GSHHS), (size_t)1, fp);
}
geo_import.end_surface();
geo_import.end_attractor();
}
// ************** MAIN PLUGIN **************
......@@ -416,10 +838,9 @@ StringXNumber GSHHSOptions_Number[] = {
{GMSH_FULLRC, "UTMZone", NULL, 0},
{GMSH_FULLRC, "UTMEquatorialRadius", NULL, 6378137},
{GMSH_FULLRC, "UTMPolarRadius", NULL, 6356752.3142},
{GMSH_FULLRC, "UTMScale", NULL, 1},
{GMSH_FULLRC, "UTMShiftX", NULL, 0},
{GMSH_FULLRC, "UTMShiftY", NULL, 0},
{GMSH_FULLRC, "WritePolarSphere",NULL,1}
{GMSH_FULLRC, "radius", NULL,6371009},
{GMSH_FULLRC, "WritePolarSphere",NULL,1},
{GMSH_FULLRC, "MinStraitsFactor",NULL,1}
};
StringXString GSHHSOptions_String[] = {
{GMSH_FULLRC, "InFileName", NULL, "gshhs_c.b"},
......@@ -442,22 +863,35 @@ void GMSH_GSHHSPlugin::getName(char *name) const
}
void GMSH_GSHHSPlugin::getInfos(char *author, char *copyright,
char *help_text) const
char *help_text) const
{
strcpy(author, "J. Lambrechts");
strcpy(copyright, "C. Geuzaine, J.-F. Remacle");
strcpy(help_text, "Plugin(GSHHS) read differenct kind of contour lines data and write a .geo file.\n"
"Valid values for \"Format\" are :\n"
" -\"gshhs\" : open GSHHS file\n"
" -\"loops2\" : import 2D contour lines in simple text format :\n"
" NB_POINTS_IN_FIRST_LOOP\n"
" COORD1 COORD2\n"
" COORD1 COORD2\n"
" ... ...\n"
" NB_POINTS_IN_SECOND_LOOP\n"
" ...\n"
);
strcpy(help_text,
"Plugin(GSHHS) read different kind of contour lines data and write a .geo file on the surface of a sphere (the Earth).\n"
"The principal application is to load GSHHS data\n (see http://www.soest.hawaii.edu/wessel/gshhs/gshhs.html).\n"
"Valid values for \"Format\" are ):\n"
" -\"gshhs\" : open GSHHS file\n"
" -\"loops2\" : import 2D contour lines in simple text format :\n"
" NB_POINTS_IN_FIRST_LOOP FIRST_LOOP_IS_CLOSED\n"
" COORD1 COORD2\n"
" COORD1 COORD2\n"
" ... ...\n"
" NB_POINTS_IN_SECOND_LOOP SECOND_LOOP_IS_CLOSED\n"
" ...\n"
" (LOOP_IS_CLOSED specify if this coast line describe a closed curve (0=no, 1=yes).\n"
"In the case of \"loops2\" format, you can specify the the coordinate system used in the input file with the"
"\"Coordinate\" option, valid values are\n"
" -\"lonlat\" for longitude-latidute radian,\n"
" -\"lonlat_degrees\" for longitude-latitude degrees,\n"
" -\"UTM\" for universal transverse mercartor (\"UTMZone\" option should be specified)\n"
" -\"cartesian\" for full 3D coordinates\n"
"\"radius\" specify the earth radius.\n"
"If the \"iField\" option is set, consecutive points closer than the value of the field iField (in meters) will not be added.\n"
"If \"MinStraitsFactor\" >0 and if a field iField is provided, coastlines closer than MinStraitsFactor*field(IField) are merged and inner corners which form an angle < pi/3 are removed.\n"
"The output is always in stereographic coordinates, if the \"WritePolarSphere\" option is not 0, a sphere is added to the geo file.\n"
"WARNING : this plugin is still experimental and need polishing and error-handling. In particular, it will probably crash if an inexistant field id is given or if the input/output cannot be open."
);
}
int GMSH_GSHHSPlugin::getNbOptions() const
......@@ -493,26 +927,26 @@ PView *GMSH_GSHHSPlugin::execute(PView * v)
char *outfilename = (char *)GSHHSOptions_String[1].def;
std::string format(GSHHSOptions_String[2].def);
std::string coordinate_name(GSHHSOptions_String[3].def);
double scale = GSHHSOptions_Number[4].def;
double shiftx = GSHHSOptions_Number[5].def;
double shifty = GSHHSOptions_Number[6].def;
bool write_polar_sphere = (bool)GSHHSOptions_Number[7].def;
int utm_zone=(int)GSHHSOptions_Number[1].def;
double utm_equatorial_radius=(double)GSHHSOptions_Number[2].def;
double utm_polar_radius=(double)GSHHSOptions_Number[3].def;
coordinate_lonlat lonlat;
coordinate_lonlat_degrees lonlat_degrees;
coordinate_utm utm(utm_zone, utm_equatorial_radius, utm_polar_radius);
GeoEarthImport geo_import(outfilename,write_polar_sphere);
double radius=(double)GSHHSOptions_Number[4].def;
bool write_polar_sphere = (bool)GSHHSOptions_Number[5].def;
double straits_factor = (double)GSHHSOptions_Number[6].def;
coordinate_lonlat lonlat(radius);
coordinate_lonlat_degrees lonlat_degrees(radius);
coordinate_utm utm(utm_zone, radius,utm_equatorial_radius, utm_polar_radius);
GeoEarthImport geo_import(outfilename,write_polar_sphere,radius);
coordinate_system *c_syst=NULL;
if(coordinate_name == "lonlat")
geo_import.set_coordinate_system(&lonlat);
c_syst=&lonlat;
else if(coordinate_name == "lonlat_degrees")
geo_import.set_coordinate_system(&lonlat_degrees);
c_syst=&lonlat_degrees;
else if(coordinate_name == "utm")
geo_import.set_coordinate_system(&utm);
c_syst=&utm;
else if(coordinate_name != "cartesian"){
Msg::Error("gshhs: Unknown coordinate system %s.\n",
coordinate_name.c_str());
coordinate_name.c_str());
return NULL;
}
Field *field = NULL;
......@@ -521,35 +955,69 @@ PView *GMSH_GSHHSPlugin::execute(PView * v)
if(!field){
Msg::Error("Field[%d] does not exist", iField);
return NULL;
}else{
geo_import.set_size_field(field);
}
}
FILE *fp;
if ((fp = fopen (filename, "rb")) == NULL ) {
Msg::Error("gshhs: Could not find file %s.\n", filename);
return NULL;
}
double x,y,z;
if(format == "gshhs") {
geo_import.set_coordinate_system(&lonlat);
import_gshhs(fp, geo_import);
SPoint3 p;
reader *read;
if(format == "loops2") {
read=new reader_loops2(filename);
}else if(format == "gshhs") {
c_syst=&lonlat;
read=new reader_gshhs(filename);
}
else if(format == "loops2") {
int npoints_in_loop;
int closed;
while(fscanf(fp, "%d %d", &npoints_in_loop,&closed) == 2) {
for(int i = 0; i < npoints_in_loop; i++) {
if(fscanf(fp, "%le %le", &x, &y) != 2) {
Msg::Error("gshhs: Error reading loops2 file \'%s\'.\n", filename);
return NULL;
}
geo_import.add_point(SPoint3(x,y,0));
loops ll;
bool closed;
while(read->next_loop(closed)!=0){
loop l;
l.closed=closed;
point oldp(0,0,0,0.);
while(read->next_point(p)){
if(c_syst)
c_syst->to_cartesian(p,p);
point newp(p[0],p[1],p[2],field,straits_factor);
if(newp.min_dist<0){
while(!l.empty() && l.back().dist(l.front())<l.back().min_dist)
l.pop_back();
l.closed=false;
if(l.size()>=3)
ll.push_back(l);
l.clear();
} else if (l.empty() || oldp.dist(newp)>newp.min_dist){
l.push_back(newp);
oldp=newp;
}
geo_import.end_loop(closed);
}
geo_import.end_surface();
geo_import.end_attractor();
}
while(!l.empty() && l.back().dist(l.front())<l.back().min_dist)
l.pop_back();
if(l.size()>=3)
ll.push_back(l);
}
delete read;
if(straits_factor>0 && iField !=0){
box *b=new box(point(-radius,-radius,-radius,0.),point(radius,radius,radius,0.));
for(loops::iterator il=ll.begin();il!=ll.end();il++)
loop_fill_box(&*il,*b);
for(loops::iterator il=ll.begin();il!=ll.end();il++)
loop_check_close_points(&*il,*b);
delete b;
for(loops::iterator il=ll.begin();il!=ll.end();il++){
box b(point(-radius,-radius,-radius,0.),point(radius,radius,radius,0.));
loop_fill_box(&*il,b);
while(false
|| loop_check_small_angles(&*il)
|| loop_check_close_points_self(&*il,b)
|| loop_check_intersections(&*il,b)
);
}
}
for(std::list<loop>::iterator l=ll.begin();l!=ll.end();l++){
for(loop::iterator p=l->begin();p!=l->end();p++)
geo_import.add_point(p->v);
geo_import.end_loop(l->closed);
}
geo_import.end_surface();
geo_import.end_attractor();
return NULL;
}
Plugin/GSHHS.h
+ 0
− 12
View file @ f6fcc770
......@@ -13,17 +13,5 @@ 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/PluginManager.cpp
+ 2
− 2
View file @ f6fcc770
......@@ -217,13 +217,13 @@ void GMSH_PluginManager::registerDefaultPlugins()
("FieldView", GMSH_RegisterFieldViewPlugin()));
allPlugins.insert(std::pair<const char*, GMSH_Plugin*>
("Triangulate", GMSH_RegisterTriangulatePlugin()));
allPlugins.insert(std::pair<const char*, GMSH_Plugin*>
("GSHHS", GMSH_RegisterGSHHSPlugin()));
#if defined(HAVE_MATH_EVAL)
allPlugins.insert(std::pair<const char*, GMSH_Plugin*>
("Evaluate", GMSH_RegisterEvaluatePlugin()));
allPlugins.insert(std::pair<const char*, GMSH_Plugin*>
("CutParametric", GMSH_RegisterCutParametricPlugin()));
allPlugins.insert(std::pair<const char*, GMSH_Plugin*>
("GSHHS", GMSH_RegisterGSHHSPlugin()));
#endif
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment