Skip to content
Snippets Groups Projects
Commit f34a3870 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

remove GSHHS plugin

parent 0bf4eaa2
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Plugin/CMakeLists.txt
+ 1
− 1
View file @ f34a3870
...@@ -15,7 +15,7 @@ set(SRC ...@@ -15,7 +15,7 @@ set(SRC
LongitudeLatitude.cpp LongitudeLatitude.cpp
Triangulate.cpp Tetrahedralize.cpp Triangulate.cpp Tetrahedralize.cpp
Warp.cpp SphericalRaise.cpp Warp.cpp SphericalRaise.cpp
Skin.cpp GSHHS.cpp Skin.cpp
MathEval.cpp ModifyComponent.cpp ExtractElements.cpp MathEval.cpp ModifyComponent.cpp ExtractElements.cpp
MakeSimplex.cpp MakeSimplex.cpp
Integrate.cpp Gradient.cpp Curl.cpp Divergence.cpp Integrate.cpp Gradient.cpp Curl.cpp Divergence.cpp
......
Plugin/GSHHS.cpp deleted 100644 → 0
+ 0
− 1095
View file @ 0bf4eaa2
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sstream>
#include <fstream>
#include "SPoint2.h"
#include "GSHHS.h"
#include "GModel.h"
#if defined(HAVE_MESH)
#include "Field.h"
#else
class Field {
public:
virtual double operator() (double x, double y, double z){ return 0.; }
};
#endif
#if defined(_MSC_VER)
inline double round(double x) { return floor(x + 0.5); }
#endif
class GMSH_GSHHSPlugin : public GMSH_PostPlugin
{
public:
// ************** Inputs (readers) *************
class reader{
public:
virtual int next_loop(bool &closed)=0;
virtual bool next_point(SPoint3 &point)=0;
virtual ~reader(){}
};
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");
}
virtual ~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.");
return true;
}
};
class reader_gshhs : public reader{
/* $Id: GSHHS.cpp,v 1.34 2009-10-12 15:46:20 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.
*/
#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
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 */
int area_full; /* Area of original full-resolution polygon in 1/10 km^2 */
int container; /* Id of container polygon that encloses this polygon (-1 if none) */
int ancestor; /* Id of ancestor polygon in the full resolution set that was the source of this polygon (-1 if none) */
};
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(),"rb");
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 (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(POINT)), SEEK_CUR);
if(first_loop)
first_loop=false;
else
max_east = 180000000; /* Only Eurasiafrica needs 270 */
}
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(POINT), (size_t)1, fp) != 1) {
Msg::Error("gshhs: Error reading gshhs file.");
return false;
}
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;
virtual ~coordinate_system(){}
};
// ************** 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);
}
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) : b(0), min_dist(0.), l(0)
{
v[0]=_x; v[1]=_y; v[2]=_z;
if(f)
min_dist=(*f)(v[0],v[1],v[2]);
}
point(double _x,double _y,double _z,double _min_dist) : b(0), l(0)
{
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)
{
double r=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
if(inverse_stereo){
xp=v[1]/(r-v[2]);
yp=v[0]/(r-v[2]);
}
else{
xp=-v[0]/(r+v[2]);
yp=-v[1]/(r+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);
}
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);
}
// 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=false)
{
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 (int)round(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]);
}
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;
}
void loop_fill_box(loop *l,box &b)
{
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;
}
bool loop_check_close_points_self(loop *l,box &b)
{
bool result=false;
int orientation = l->orientation(l->begin(), --l->end());
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)!=orientation || l->length(i,id1)<3)){
i=l->remove_range(i,id1);
result=true;
}
if(not_a_loop!=1 && (l->orientation(id1,i)!=orientation || l->length(id1,i)<3)){
i=l->remove_range(id1,i);
result=true;
}else
i++;
}else
i++;
}
return result;
}
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;
}
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];
}
l->insert_range(++f0,f1->l->next(f1),f1);
cp[j]->l->clear();
merged=true;
break;
}
}
if(!merged)
i++;
}
return true;
}
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;
}
int _write_polar_sphere;
public:
GeoEarthImport(const std::string _filename, int write_polar_sphere,double radius)
{
_write_polar_sphere = write_polar_sphere;
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 > 0){
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)
{
double r=sqrt(point.x()*point.x()+point.y()*point.y()+point.z()*point.z());
SPoint2 stereo(-point.x() / (r + point.z()), -point.y() / (r + point.z()));
if (_write_polar_sphere == -2){
stereo = SPoint2(2 * r * point.x() / (r + point.z()), 2 * r * point.y() / (r + point.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 << "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;
}
}
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();
}
};
std::string getName() const { return "GSHHS"; }
std::string getShortHelp() const
{
return "Import and process GSHHS data sets";
}
std::string getHelp() const;
std::string getAuthor() const { return "J. Lambrechts"; }
int getNbOptions() const;
int getNbOptionsStr() const;
StringXNumber *getOption(int iopt);
StringXString *getOptionStr(int iopt);
PView *execute(PView *);
};
// ************** MAIN PLUGIN **************
StringXNumber GSHHSOptions_Number[] = {
{GMSH_FULLRC, "iField", NULL, -1.},
{GMSH_FULLRC, "UTMZone", NULL, 0},
{GMSH_FULLRC, "UTMEquatorialRadius", NULL, 6378137},
{GMSH_FULLRC, "UTMPolarRadius", NULL, 6356752.3142},
{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"},
{GMSH_FULLRC, "OutFileName", NULL, "earth.geo"},
{GMSH_FULLRC, "Format", NULL, "gshhs"},
{GMSH_FULLRC, "Coordinate", NULL, "cartesian"}
};
extern "C"
{
GMSH_Plugin *GMSH_RegisterGSHHSPlugin()
{
return new GMSH_GSHHSPlugin();
}
}
std::string GMSH_GSHHSPlugin::getHelp() const
{
return
"Plugin(GSHHS) read different kind of contour lines data "
"and write a .geo file on the surface of a sphere (the Earth).\n\n"
"The principal application is to load GSHHS data\n (see "
"http://www.soest.hawaii.edu/wessel/gshhs/gshhs.html).\n\n"
"Valid values for \"Format\" are:\n\n"
"- \"gshhs\": open GSHHS file\n\n"
"- \"loops2\": import 2D contour lines in simple text format:\n\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\n"
"(LOOP_IS_CLOSED specifies if this coast line describes a closed "
"curve (0=no, 1=yes)).\n\n"
"In the case of \"loops2\" format, you "
"can specify the coordinate system used in the input file "
"with the \"Coordinate\" option. Valid values are\n\n"
"- \"lonlat\" for longitude-latidute radian,\n\n"
"- \"lonlat_degrees\" for longitude-latitude degrees,\n\n"
"- \"UTM\" for universal transverse mercartor (\"UTMZone\" "
"option should be specified)\n\n"
"- \"cartesian\" for full 3D coordinates\n\n"
"- \"radius\" specify the earth radius.\n\n"
"If the \"iField\" option is set, consecutive points closer "
"than the value of the field iField (in meters) will not be "
"added.\n\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\n"
"The output is always in stereographic coordinates, if "
"the \"WritePolarSphere\" option is greater than 0, a sphere is "
"added to the geo file.\n\n"
"WARNING: this plugin is still experimental and needs "
"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
{
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];
}
PView *GMSH_GSHHSPlugin::execute(PView * v)
{
void projector(SPoint2,SPoint3);
int iField = (int)GSHHSOptions_Number[0].def;
char *filename = (char *)GSHHSOptions_String[0].def.c_str();
char *outfilename = (char *)GSHHSOptions_String[1].def.c_str();
std::string format(GSHHSOptions_String[2].def);
std::string coordinate_name(GSHHSOptions_String[3].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;
double radius=(double)GSHHSOptions_Number[4].def;
int write_polar_sphere = (int)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")
c_syst=&lonlat;
else if(coordinate_name == "lonlat_degrees")
c_syst=&lonlat_degrees;
else if(coordinate_name == "utm")
c_syst=&utm;
else if(coordinate_name != "cartesian"){
Msg::Error("gshhs: Unknown coordinate system %s.\n",
coordinate_name.c_str());
return NULL;
}
Field *field = NULL;
#if defined(HAVE_MESH)
if (iField != -1) {
field = GModel::current()->getFields()->get(iField);
if(!field){
Msg::Error("Field[%d] does not exist", iField);
return NULL;
}
}
#endif
SPoint3 p;
reader *read=0;
if(format == "loops2") {
read=new reader_loops2(filename);
}
else if(format == "gshhs") {
c_syst=&lonlat;
read=new reader_gshhs(filename);
}
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);
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;
}
}
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){
for(loops::iterator il=ll.begin();il!=ll.end();il++)
for(loop::iterator ip=il->begin();ip!=il->end();ip++)
ip->min_dist*=straits_factor;
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 deleted 100644 → 0
+ 0
− 16
View file @ 0bf4eaa2
// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
#ifndef _GSHHS_H_
#define _GSHHS_H_
#include "Plugin.h"
extern "C"
{
GMSH_Plugin *GMSH_RegisterGSHHSPlugin();
}
#endif
Plugin/PluginManager.cpp
+ 0
− 3
View file @ f34a3870
...@@ -47,7 +47,6 @@ ...@@ -47,7 +47,6 @@
#include "Lambda2.h" #include "Lambda2.h"
#include "ModifyComponent.h" #include "ModifyComponent.h"
#include "Probe.h" #include "Probe.h"
#include "GSHHS.h"
#include "HomologyComputation.h" #include "HomologyComputation.h"
#include "HomologyPostProcessing.h" #include "HomologyPostProcessing.h"
#include "ExtractEdges.h" #include "ExtractEdges.h"
...@@ -227,8 +226,6 @@ void PluginManager::registerDefaultPlugins() ...@@ -227,8 +226,6 @@ void PluginManager::registerDefaultPlugins()
("Probe", GMSH_RegisterProbePlugin())); ("Probe", GMSH_RegisterProbePlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*> allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("Triangulate", GMSH_RegisterTriangulatePlugin())); ("Triangulate", GMSH_RegisterTriangulatePlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("GSHHS", GMSH_RegisterGSHHSPlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*> allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("ExtractEdges", GMSH_RegisterExtractEdgesPlugin())); ("ExtractEdges", GMSH_RegisterExtractEdgesPlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*> allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment