## @namespace Geo::OGC
# @brief The simple feature geometries
#
# The classes and methods in the Geo::OGC:: namespace should conform
# to the OGC implementation specification (currently 06-103r3) for
# Geographic Information - Simple feature access.
# @note Some of the brief descriptions in this document are taken
# verbatim from <a
# href="http://www.opengeospatial.org/standards/sfa">the OGC
# specification (OGC 06-103r3 Version 1.2.0)</a>. <a
# href="http://www.opengeospatial.org/ogc/document">Copyright (c) 2006
# Open Geospatial Consortium, Inc. All Rights Reserved.</a>
# @note Calling a method that is not yet implemented causes an
# exception. The status of implementation is not always shown in this
# document.
# @note Most of the methods for testing spatial relations and the
# methods that support spatial analysis are not yet implemented.
package Geo::OGC::Geometry;
=pod
=head1 NAME
Geo::OGC::Geometry - Simple feature geometry classes
The classes and methods in the Geo::OGC:: namespace should conform to
the OGC (opengeospatial.org) implementation specification (currently
06-103r3) for Geographic Information - Simple feature access.
This module is documented using the doxygen format. Documentation in
HTML and in other formats can be generated with <a
href="http://www.stack.nl/~dimitri/doxygen/">doxygen</a> and <a
href="http://www.bigsister.ch/doxygenfilter/">perl doxygen filter</a>.
The latest version of the documentation is automatically generated at
http://map.hut.fi/doc/Geoinformatica/html/class_geo_1_1_o_g_c_1_1_geometry.html
=cut
use strict;
use Carp;
# Force to "C" locale
use POSIX;
POSIX::setlocale( &POSIX::LC_NUMERIC, "C" );
BEGIN {
use Exporter 'import';
use vars qw /$SNAP_DISTANCE_SQR/;
our %EXPORT_TAGS = ( 'all' => [ qw( &ccw &intersect &distance_point_line_sqr
$SNAP_DISTANCE_SQR ) ] );
our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );
our $VERSION = '0.05';
$SNAP_DISTANCE_SQR = 1E-6;
}
## @cmethod Geo::OGC::Geometry new(%params)
# @brief Create a new and initialized object.
#
# @param params A set of named parameters (the subclasses may add new
# known parameters):
# - <i>Text</i> A well-known text, constructs an object of class the
# text defines.
# - <i>SRID</i> Sets the SRID attribute of the object, default is -1.
# - <i>Precision</i> If specified, has the effect that numeric
# comparisons in the Equals method is is preceded with a rounding
# operation (using sprintf "%.pe", where p is the Precision-1, i.e the
# number of meaningful numbers is Precision). Affects also AsText.
#
# This method should be called eventually by all constructors. Blesses
# the object into the correct class and calls init.
sub new {
my($package, %params) = @_;
return parse_wkt($params{Text}) if exists $params{Text};
my $self = {};
bless $self => (ref($package) or $package);
$self->init(%params);
return $self;
}
## @method void init(%params)
# @brief Override in new classes but call $self->SUPER::init(%params);
sub init {
my($self, %params) = @_;
$self->{SRID} = exists $params{SRID} ? $params{SRID} : -1;
$self->{Precision} = $params{Precision}-1 if exists $params{Precision};
}
## @method void copy($clone)
# @brief Copies the attributes of this object to the other object,
# which is a newly created object of same class.
sub copy {
my($self, $clone) = @_;
$clone->{SRID} = $self->{SRID};
$clone->{Precision} = $self->{Precision} if exists $self->{Precision};
}
## @method Geo::OGC::Geometry parse_wkt($Text)
# @brief parse well known text and construct respective geometry
sub parse_wkt {
my $self;
#print STDERR "parse: $_[0]\n";
for ($_[0]) {
if (/^\s*POINT/i) {
s/^\s*POINT\s*([ZM]*)\s*\(\s*//i;
my $m = lc($1);
s/\s*\)\s*$//;
my @point = split /\s+/;
$self = Geo::OGC::Point->new("point$m"=>\@point);
} elsif (/^\s*MULTIPOINT/i) {
s/^\s*MULTIPOINT\s*([ZM]*)\s*\(\s*//i;
my $m = $1;
s/\s*\)\s*$//;
my @points = split /\s*,\s*/;
$self = Geo::OGC::MultiPoint->new();
for my $p (@points) {
$self->AddGeometry(parse_wkt("POINT$m ($p)"));
}
} elsif (/^\s*LINESTRING/i) {
s/^\s*LINESTRING\s*([ZM]*)\s*\(\s*//i;
my $m = $1;
s/\s*\)\s*$//;
my @points = split /\s*,\s*/;
$self = Geo::OGC::LineString->new();
for my $p (@points) {
$self->AddPoint(parse_wkt("POINT$m ($p)"));
}
} elsif (/^\s*MULTILINESTRING/i) {
s/^\s*MULTILINESTRING\s*([ZM]*)[\s\(]*//i;
my $m = $1;
s/[\s\)]*$//;
my @strings = split /\)\s*,\s*\(/;
$self = Geo::OGC::MultiLineString->new();
for my $s (@strings) {
$self->AddGeometry(parse_wkt("LINESTRING$m ($s)"));
}
} elsif (/^\s*LINEARRING/i) {
s/^\s*LINEARRING\s*([ZM]*)\s*\(\s*//i;
my $m = $1;
s/\s*\)\s*$//;
my @points = split /\s*,\s*/;
$self = Geo::OGC::LinearRing->new();
for my $p (@points) {
$self->AddPoint(parse_wkt("POINT$m ($p)"));
}
} elsif (/^\s*POLYGON/i) {
s/^\s*POLYGON\s*([ZM]*)[\s\(]*//i;
my $m = $1;
s/[\s\)]*$//;
my @rings = split /\)\s*,\s*\(/;
$self = Geo::OGC::Polygon->new();
$self->ExteriorRing(parse_wkt("LINEARRING$m (".shift(@rings).")"));
for my $r (@rings) {
$self->AddInteriorRing(parse_wkt("LINEARRING$m ($r)"));
}
} elsif (/^\s*POLYHEDRALSURFACE/i) {
s/^\s*POLYHEDRALSURFACE\s*([ZM]*)[\s\(]*//i;
my $m = $1;
s/[\s\)]*$//;
my @patches = split /\)\s*,\s*\(/;
$self = Geo::OGC::PolyhedralSurface->new();
for my $p (@patches) {
$self->AddPatch(parse_wkt("POLYGON$m (($p)"));
}
} elsif (/^\s*MULTIPOLYGON/i) {
s/^\s*MULTIPOLYGON\s*([ZM]*)[\s\(]*//i;
my $m = $1;
s/[\s\)]*$//;
my @polygons = split /\)\s*\)\s*,\s*\(\s*\(/;
$self = Geo::OGC::MultiPolygon->new();
for my $p (@polygons) {
$self->AddGeometry(parse_wkt("POLYGON$m (($p))"));
}
} elsif (/^\s*GEOMETRYCOLLECTION/i) {
s/^\s*GEOMETRYCOLLECTION\s*([ZM]*)\s*\(\s*//i;
my $m = $1;
s/\s*\)\s*$//;
my @b = /,([A-Z])/g;
unshift @b,'';
my @geometries = split /,[A-Z]/;
$self = Geo::OGC::GeometryCollection->new();
for my $i (0..$#geometries) {
$self->AddGeometry(parse_wkt($b[$i].$geometries[$i]));
}
} else {
my $beginning = substr $_, 0, 20;
croak "can't parse text beginning as: $beginning";
}
}
return $self;
}
## @fn $ccw($x0, $y0, $x1, $y1, $x2, $y2)
# @brief counterclockwise from Sedgewick: Algorithms in C
# @return -1, 0, or 1
sub ccw {
my($x0, $y0, $x1, $y1, $x2, $y2) = @_;
my $dx1 = $x1 - $x0;
my $dy1 = $y1 - $y0;
my $dx2 = $x2 - $x0;
my $dy2 = $y2 - $y0;
return +1 if $dx1*$dy2 > $dy1*$dx2;
return -1 if $dx1*$dy2 < $dy1*$dx2;
return -1 if ($dx1*$dx2 < 0) or ($dy1*$dy2 < 0);
return +1 if ($dx1*$dx1+$dy1*$dy1) < ($dx2*$dx2+$dy2*$dy2);
return 0;
}
## @fn $intersect($x11, $y11, $x12, $y12, $x21, $y21, $x22, $y22)
# @brief Test intersection of two lines from Sedgewick: Algorithms in C
sub intersect {
my($x11, $y11, $x12, $y12, $x21, $y21, $x22, $y22) = @_;
return ((ccw($x11, $y11, $x12, $y12, $x21, $y21)
*ccw($x11, $y11, $x12, $y12, $x22, $y22)) <= 0)
&& ((ccw($x21, $y21, $x22, $y22, $x11, $y11)
*ccw($x21, $y21, $x22, $y22, $x12, $y12)) <= 0);
}
sub intersection_point {
my($x11, $y11, $x12, $y12, $x21, $y21, $x22, $y22) = @_;
my $dy1 = $y12 - $y11;
my $dx1 = $x12 - $x11;
my $dy2 = $y22 - $y21;
my $dx2 = $x22 - $x21;
# (dy1*dx2 - dy2*dx1)*x = dx1*dx2*(y21-y11) - dy2*dx1*x21 + dy1*dx2*x11
# (dy1*dx1 - dy2*dx2)*y = dy1*dy2*(x21-x11) - dy1*dx2*y21 + dy2*dx1*y11
my $x = ($dx1*$dx2*($y21-$y11) - $dy2*$dx1*$x21 + $dy1*$dx2*$x11)/($dy1*$dx2 - $dy2*$dx1);
my $y = ($dy1*$dy2*($x21-$x11) - $dy1*$dx2*$y21 + $dy2*$dx1*$y11)/($dy1*$dx1 - $dy2*$dx2);
return ($x, $y);
}
## @fn $distance_point_line_sqr($x, $y, $x1, $y1, $x2, $y2)
# @brief Compute the distance of a point to a line.
sub distance_point_line_sqr {
my($x, $y, $x1, $y1, $x2, $y2) = @_;
my $dx = $x2-$x1;
my $dy = $y2-$y1;
my $l2 = $dx*$dx + $dy*$dy;
my $u = (($x - $x1) * $dx + ($y - $y1) * $dy) / $l2;
if ($u < 0) { # distance to point 1
return (($x-$x1)*($x-$x1) + ($y-$y1)*($y-$y1));
} elsif ($u > 1) { # distance to point 2
return (($x-$x2)*($x-$x2) + ($y-$y2)*($y-$y2));
} else {
my $ix = $x1 + $u * $dx;
my $iy = $y1 + $u * $dy;
return (($x-$ix)*($x-$ix) + ($y-$iy)*($y-$iy));
}
}
## @method Geo::OGC::Geometry Clone()
# @brief Clones this object, i.e., creates a spatially exact copy.
sub Clone {
my($self) = @_;
my $clone = Geo::OGC::Geometry::new($self);
$self->copy($clone);
return $clone;
}
## @method $Precision($Precision)
# @brief Sets or gets the precision
# @note Not in the specification
sub Precision {
my($self, $Precision) = @_;
defined $Precision ?
$self->{Precision} = $Precision-1 : $self->{Precision}+1;
}
## @method $Dimension()
# @brief The dimension of this geometric object. In non-homogeneous
# collections, this will return the largest topological dimension of
# the contained objects.
# @return 2 or 3
sub Dimension {
my($self) = @_;
croak "Dimension method for class ".ref($self)." is not implemented yet";
}
## @method $GeometryType()
# @brief Returns the name of the class of this geometric object.
sub GeometryType {
my($self) = @_;
croak "GeometryType method for class ".ref($self)." is not implemented yet";
}
## @method $SRID($SRID)
# @brief Returns (or sets) the Spatial Reference System ID for this
# geometric object.
# @param SRID [optional] The new SRID if setting it.
sub SRID {
my($self, $SRID) = @_;
defined $SRID ?
$self->{SRID} = $SRID : $self->{SRID};
}
## @method Geo::OGC::Polygon Envelope()
# @brief The minimum bounding box for this Geometry.
# @note This library returns always the envelope as a ring [(minx,
# miny), (maxx, miny), (maxx, maxy), (minx, maxy), (minx, miny)]
sub Envelope {
my($self) = @_;
croak "Envelope method for class ".ref($self)." is not implemented yet";
}
## @method $as_text($force_parens, $include_tag)
# @brief A helper method used by AsText
sub as_text {
my($self, $force_parens, $include_tag) = @_;
croak "as_text method for class ".ref($self)." is not implemented yet";
}
## @method $AsText()
# @brief Returns this object in Well-known Text Representation of Geometry
sub AsText {
my($self) = @_;
return $self->as_text(1, 1);
}
## @method $AsBinary()
# @brief Returns this object in Well-known Binary Representation of Geometry
# @note Not implemented yet.
sub AsBinary {
my($self) = @_;
croak "AsBinary method for class ".ref($self)." is not implemented yet";
}
## @method $IsEmpty()
# @brief Returns true if this object is empty, i.e. contains no points
# or no data.
sub IsEmpty {
my($self) = @_;
croak "IsEmpty method for class ".ref($self)." is not implemented yet";
}
## @method $IsSimple()
# @brief Returns true if this geometric object has no anomalous
# geometric points, such as self intersection or self tangency.
sub IsSimple {
my($self) = @_;
croak "IsSimple method for class ".ref($self)." is not implemented yet";
}
## @method $Is3D()
# @brief Returns true if this geometric object has z coordinate values.
sub Is3D {
my($self) = @_;
croak "Is3D method for class ".ref($self)." is not implemented yet";
}
## @method $IsMeasured()
# @brief Returns true if this geometric object has m coordinate
# values.
sub IsMeasured {
my($self) = @_;
croak "IsMeasured method for class ".ref($self)." is not implemented yet";
}
## @method Geo::OGC::Geometry Boundary()
# @brief Returns the closure of the combinatorial boundary of this
# geometric object.
# @note Not implemented yet.
sub Boundary {
my($self) = @_;
croak "Boundary method for class ".ref($self)." is not implemented yet";
}
## @method $Equals($geom)
# @brief Returns true if this geometric object is "spatially equal" to
# the another geometry.
sub Equals {
my($self, $geom) = @_;
croak "Equals method for class ".ref($self)." is not implemented yet";
}
sub Disjoint {
my($self, $geom) = @_;
croak "Disjoint method for class ".ref($self)." is not implemented yet";
}
sub Intersects {
my($self, $geom) = @_;
croak "Intersects method for class ".ref($self)." is not implemented yet";
}
sub Touches {
my($self, $geom) = @_;
croak "Touches method for class ".ref($self)." is not implemented yet";
}
sub Crosses {
my($self, $geom) = @_;
croak "Crosses method for class ".ref($self)." is not implemented yet";
}
sub Within {
my($self, $geom) = @_;
croak "Within method for class ".ref($self)." is not implemented yet";
}
sub Contains {
my($self, $geom) = @_;
croak "Contains method for class ".ref($self)." is not implemented yet";
}
sub Overlaps {
my($self, $geom) = @_;
croak "Overlaps method for class ".ref($self)." is not implemented yet";
}
sub Relate {
my($self, $geom, $int_pat) = @_;
croak "Relate method for class ".ref($self)." is not implemented yet";
}
sub LocateAlong {
my($self, $mValue) = @_;
croak "LocateAlong method for class ".ref($self)." is not implemented yet";
}
sub LocateBetween {
my($self, $mStart, $mEnd) = @_;
croak "LocateBetween method for class ".ref($self)." is not implemented yet";
}
sub Distance {
my($self, $geom) = @_;
croak "Distance method for class ".ref($self)." is not implemented yet";
}
sub Buffer {
my($self, $distance) = @_;
croak "Buffer method for class ".ref($self)." is not implemented yet";
}
sub ConvexHull {
my($self) = @_;
croak "ConvexHull method for class ".ref($self)." is not implemented yet";
}
sub Intersection {
my($self, $geom) = @_;
croak "Intersection method for class ".ref($self)." is not implemented yet";
}
sub Union {
my($self, $geom) = @_;
croak "Union method for class ".ref($self)." is not implemented yet";
}
sub Difference {
my($self, $geom) = @_;
croak "Difference method for class ".ref($self)." is not implemented yet";
}
sub SymDifference {
my($self, $geom) = @_;
croak "SymDifference method for class ".ref($self)." is not implemented yet";
}
## @method MakeCollection()
# @brief Creates a collection which contains this geometry
# @note Not in the specification
sub MakeCollection {
my($self) = @_;
croak "MakeCollection method for class ".ref($self)." is not implemented";
}
## @method ApplyTransformation($transf)
# @param transf A point transformation method which will be applied
# for all points in the geometry as:
# @code
# ($new_x, $new_y, $new_z) = $transf->($x, $y, $z)
# @endcode
# @note Not in the specification
sub ApplyTransformation {
my($self, $transf) = @_;
croak "ApplyTransformation method for class ".ref($self)." is not implemented";
}
## @method LastPolygon()
# @brief Returns last (latest added) polygon or undef
sub LastPolygon {
return undef;
}
#
# SpatialReferenceSystem
#
package Geo::OGC::SpatialReferenceSystem;
use strict;
use Carp;
sub new {
my($package, %params) = @_;
my $self = {};
bless $self => (ref($package) or $package);
}
#
# Point
#
package Geo::OGC::Point;
use strict;
use Carp;
use Geo::OGC::Geometry qw/:all/;
our @ISA = qw( Geo::OGC::Geometry );
## @cmethod new(%params)
# @brief Construct a new point
# @param params The following syntaxes are allowed:
# @code
# $point = Geo::OGC::Point->new($x, $y);
# $point = Geo::OGC::Point->new($x, $y, $z);
# $point = Geo::OGC::Point->new(point => [$x, $y]);
# $point = Geo::OGC::Point->new(point => [$x, $y, $z]);
# $point = Geo::OGC::Point->new(point => [$x, $y, $z, $m]);
# $point = Geo::OGC::Point->new(pointz => [$x, $y, $z]);
# $point = Geo::OGC::Point->new(pointz => [$x, $y, $z, $m]);
# $point = Geo::OGC::Point->new(pointm => [$x, $y, $m]);
# $point = Geo::OGC::Point->new(pointm => [$x, $y, $z, $m]);
# $point = Geo::OGC::Point->new(pointzm => [$x, $y, $z, $m]);
# $point = Geo::OGC::Point->new(X => $x, Y => $y);
# $point = Geo::OGC::Point->new(X => $x, Y => $y, Z => $z);
# $point = Geo::OGC::Point->new(X => $x, Y => $y, Z => $z, M => $m);
# @endcode
sub new {
my $package = shift;
my %params;
if (@_ == 2 and !($_[0] =~ /^[XYZMpP]/)) { # allow syntax Point->new($x, $y);
$params{X} = shift;
$params{Y} = shift;
} elsif (@_ == 3) { # allow syntax Point->new($x, $y, $z);
$params{X} = shift;
$params{Y} = shift;
$params{Z} = shift;
} else {
%params = @_;
}
# support comma as decimal point, and space in numbers
for my $k (keys %params) {
if (ref($params{$k})) {
for my $p (@{$params{$k}}) {
$p =~ s/,/./g;
$p =~ s/ //g;
#print STDERR "point: $_\n";
}
} else {
$params{$k} =~ s/,/./g;
$params{$k} =~ s/ //g;
#print STDERR "point: $_ => $params{$_}\n";
}
}
my $self = Geo::OGC::Geometry::new($package, %params);
return $self;
}
sub init {
my($self, %params) = @_;
$self->SUPER::init(%params);
# +0 catches non-numeric error, if warnings are on
if ($params{point}) {
$self->{X} = $params{point}[0]+0;
$self->{Y} = $params{point}[1]+0;
$self->{Z} = $params{point}[2]+0 if @{$params{point}} > 2;
$self->{M} = $params{point}[3]+0 if @{$params{point}} > 3;
} elsif ($params{pointz}) {
$self->{X} = $params{pointz}[0]+0;
$self->{Y} = $params{pointz}[1]+0;
$self->{Z} = $params{pointz}[2]+0;
if (@{$params{pointz}} == 4) {
$self->{M} = $params{pointz}[3]+0;
}
} elsif ($params{pointm}) {
$self->{X} = $params{pointm}[0]+0;
$self->{Y} = $params{pointm}[1]+0;
if (@{$params{pointm}} == 3) {
$self->{M} = $params{pointm}[2]+0;
} elsif (@{$params{pointm}} == 4) {
$self->{Z} = $params{pointm}[2]+0;
$self->{M} = $params{pointm}[3]+0;
}
} elsif ($params{pointzm}) {
$self->{X} = $params{pointm}[0]+0;
$self->{Y} = $params{pointm}[1]+0;
$self->{Z} = $params{pointm}[2]+0;
$self->{M} = $params{pointm}[3]+0;
} else {
$self->{X} = $params{X}+0 if exists $params{X};
$self->{Y} = $params{Y}+0 if exists $params{Y};
$self->{Z} = $params{Z}+0 if exists $params{Z};
$self->{M} = $params{M}+0 if exists $params{M};
}
}
sub copy {
my($self, $clone) = @_;
$self->SUPER::copy($clone);
for my $a (qw/X Y Z M/) {
$clone->{$a} = $self->{$a} if exists $self->{$a};
}
}
## @method point()
# @brief Return a reference to an anonymous array that contains the point data.
# @note Note that there is no difference between [x,y,z] and [x,y,m]
sub point {
my($self) = @_;
my @point = ($self->{X}, $self->{Y});
push @point, $self->{Z} if exists $self->{Z};
push @point, $self->{M} if exists $self->{M};
return [@point];
}
sub GeometryType {
return 'Point';
}
sub Dimension {
return 0;
}
sub Clone {
my($self) = @_;
my $m = exists $self->{M} ? 'm' : '';
return Geo::OGC::Point::new($self, "point$m" => $self->point);
}
sub IsEmpty {
my($self) = @_;
return !(exists $self->{X});
}
## @method IsSimple()
# @brief A point is always simple.
sub IsSimple {
my($self) = @_;
return 1;
}
sub Is3D {
my($self) = @_;
return exists $self->{Z};
}
sub IsMeasured {
my($self) = @_;
return exists $self->{M};
}
sub Boundary {
my($self) = @_;
return $self->Clone;
}
sub X {
my($self, $X) = @_;
defined $X ?
$self->{X} = $X : $self->{X};
}
sub Y {
my($self, $Y) = @_;
defined $Y ?
$self->{Y} = $Y : $self->{Y};
}
## @method Z($Z)
# @brief sets or gets the z coordinate
# @param Z [optional]
# @note setting is not in the specification
# @note returns undef if z does not exist or if it exists but is undefined
sub Z {
my($self, $Z) = @_;
defined $Z ?
$self->{Z} = $Z : (exists $self->{Z} ? $self->{Z} : undef);
}
## @method M($M)
# @brief sets or gets the measure
# @param M [optional]
# @note setting is not in the specification
# @note returns undef if M does not exist or if it exists but is undefined
sub M {
my($self, $M) = @_;
defined $M ?
$self->{M} = $M : (exists $self->{M} ? $self->{M} : undef);
}
sub as_text {
my($self, $force_parens, $include_tag) = @_;
my @coords;
my $ZM = exists $self->{Z} ? 'Z' : '';
if (exists $self->{Precision}) {
for my $a (qw/X Y Z/) {
last unless exists $self->{$a};
my $s = sprintf('%.'.$self->{Precision}.'e', $self->{$a});
push @coords, $s;
}
} else {
for my $a (qw/X Y Z/) {
push @coords, $self->{$a} if exists $self->{$a};
}
}
if (exists $self->{M}) {
push @coords, $self->{M};
$ZM .= 'M';
}
my $text = join(' ', @coords);
$text =~ s/,/./g; # setting POSIX numeric locale does not seem to work??
$text = '('.$text.')' if $force_parens;
$text = "POINT$ZM ".$text if $include_tag;
return $text;
}
# what should we do with z?
sub Equals {
my($self, $geom) = @_;
return 0 unless $geom->isa('Geo::OGC::Point');
if (exists $self->{Precision}) {
for my $a (qw/X Y Z/) {
last unless exists $self->{$a} and exists $geom->{$a};
my $s = sprintf('%.'.$self->{Precision}.'e', $self->{$a});
my $g = sprintf('%.'.$self->{Precision}.'e', $geom->{$a});
return 0 if $s != $g;
}
return 1;
}
return (($self->{X}-$geom->{X})**2+($self->{Y}-$geom->{Y})**2) < $SNAP_DISTANCE_SQR;
}
sub DistanceToLineStringSqr {
my($self, $linestring) = @_;
my($x, $y) = ($self->{X}, $self->{Y});
my $p1 = $linestring->{Points}[0];
return unless $p1;
my $p2 = $linestring->{Points}[1];
return (($x-$p1->{X})**2+($y-$p1->{Y})**2) unless $p2;
my $distance = distance_point_line_sqr($x, $y, $p1->{X}, $p1->{Y}, $p2->{X}, $p2->{Y});
for my $i (2..$#{$linestring->{Points}}) {
$p1 = $p2;
$p2 = $linestring->{Points}[$i];
my $d = distance_point_line_sqr($x, $y, $p1->{X}, $p1->{Y}, $p2->{X}, $p2->{Y});
$distance = $d if $d < $distance;
}
return $distance;
}
sub Distance {
my($self, $geom) = @_;
if ($geom->isa('Geo::OGC::Point')) {
return sqrt(($self->{X}-$geom->{X})**2 + ($self->{Y}-$geom->{Y})**2);
} elsif ($geom->isa('Geo::OGC::LineString')) {
return sqrt($self->DistanceToLineStringSqr($geom));
} elsif ($geom->isa('Geo::OGC::Polygon')) {
if ($geom->{ExteriorRing}->IsPointIn($self)) {
for my $ring (@{$geom->{InteriorRings}}) {
return sqrt($self->DistanceToLineStringSqr($ring)) if $ring->IsPointIn($self);
}
return 0;
} else {
return sqrt($self->DistanceToLineStringSqr($geom->{ExteriorRing}));
}
} elsif ($geom->isa('Geo::OGC::GeometryCollection')) {
my $dist = $self->Distance($geom->{Geometries}[0]);
for my $g (@{$geom->{Geometries}}[1..$#{$geom->{Geometries}}]) {
my $d = $self->Distance($g);
$dist = $d if $d < $dist;
}
return $dist;
} else {
croak "can't compute distance between a ".ref($geom)." and a point";
}
}
# what should this be?
sub Envelope {
my($self) = @_;
my $r = Geo::OGC::LinearRing->new;
$r->AddPoint(Geo::OGC::Point->new($self->{X}, $self->{Y}));
$r->AddPoint(Geo::OGC::Point->new($self->{X}, $self->{Y}));
$r->AddPoint(Geo::OGC::Point->new($self->{X}, $self->{Y}));
$r->AddPoint(Geo::OGC::Point->new($self->{X}, $self->{Y}));
$r->Close;
return $r;
}
sub Area {
return 0;
}
sub Intersection {
my($self, $geom) = @_;
return $self->Clone if $self->Within($geom);
return undef;
}
sub Within {
my($self, $geom) = @_;
if ($geom->isa('Geo::OGC::Point')) {
return $self->Equals($geom) ? 1 : 0;
} elsif ($geom->isa('Geo::OGC::LineString')) {
return $self->DistanceToLineStringSqr($geom) < $SNAP_DISTANCE_SQR ? 1 : 0;
} elsif ($geom->isa('Geo::OGC::Polygon')) {
if (!($geom->{ExteriorRing}->IsPointStricktlyOut($self))) {
for my $ring (@{$geom->{InteriorRing}}) {
return 0 if $ring->IsPointStricktlyIn($self);
}
return 1;
}
return 0;
} elsif ($geom->isa('Geo::OGC::GeometryCollection')) {
for my $g (@{$geom->{Geometries}}) {
return 1 if $self->Within($g);
}
return 0;
} else {
croak "point within ".ref($geom)." is not implemented yet";
}
}
sub MakeCollection {
my($self) = @_;
my $coll = Geo::OGC::MultiPoint->new;
$coll->AddGeometry($self);
return $coll;
}
sub ApplyTransformation {
my($self, $transf) = @_;
if (@_ > 2) {
($self->{X}, $self->{Y}, $self->{Z}) = $transf->($self->{X}, $self->{Y}, $self->{Z});
} else {
($self->{X}, $self->{Y}) = $transf->($self->{X}, $self->{Y});
}
}
sub ClosestVertex {
my($self, $x, $y) = @_;
return (($self->{X}-$x)**2 + ($self->{Y}-$y)**2);
}
sub VertexAt {
my $self = shift;
return ($self);
}
sub ClosestPoint {
my($self, $x, $y) = @_;
return (($self->{X}-$x)**2 + ($self->{Y}-$y)**2);
}
sub AddVertex {
}
sub DeleteVertex {
}
#
# Curve
#
package Geo::OGC::Curve;
# @brief 1-dimensional geometric object
#
# Curve is implemented as a sequence of Points.
use strict;
use Carp;
use Geo::OGC::Geometry qw/:all/;
our @ISA = qw( Geo::OGC::Geometry );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::Geometry::new($package, %params);
return $self;
}
sub init {
my($self, %params) = @_;
$self->SUPER::init(%params);
$self->{Points} = [];
if ($params{points}) {
for my $p (@{$params{points}}) {
$self->AddPoint(Geo::OGC::Point->new(point=>$p));
}
} elsif ($params{pointsm}) {
for my $p (@{$params{pointsm}}) {
$self->AddPoint(Geo::OGC::Point->new(pointm=>$p));
}
}
}
sub copy {
my($self, $clone) = @_;
$self->SUPER::copy($clone);
for my $p (@{$self->{Points}}) {
$clone->AddPoint($p->Clone);
}
}
sub GeometryType {
return 'Curve';
}
sub Dimension {
return 1;
}
sub as_text {
my($self, $force_parens, $include_tag) = @_;
my $text = join(',', map {$_->as_text} @{$self->{Points}});
$text = '('.$text.')';
my $ZM = $self->Is3D ? 'Z' : '';
$ZM .= 'M' if $self->IsMeasured;
$text = uc($self->GeometryType).$ZM.' '.$text if $include_tag;
return $text;
}
## @method AddPoint($point, $i)
# @param point A Point object
# @param i [optional] The location in the sequence (1..N+1) where to add the Point.
# Adds to the end (N+1) by default.
# @note not in the specification
sub AddPoint {
my($self, $point, $i) = @_;
croak 'usage: Curve->AddPoint($point) '
unless $point and $point->isa('Geo::OGC::Point');
my $points = $self->{Points};
if (defined $i) {
my $temp = $points->[$i-1];
splice @$points, $i-1, 1, ($point, $temp);
} else {
push @$points, $point;
}
}
## @method DeletePoint($i)
# @param i The location in the sequence (1..N) from where to delete the Point.
# @note not in the specification
sub DeletePoint {
my($self, $i) = @_;
my $points = $self->{Points};
splice @$points, $i-1, 1;
}
sub StartPoint {
my($self) = @_;
my $points = $self->{Points};
return $points->[0] if @$points;
}
sub EndPoint {
my($self) = @_;
my $points = $self->{Points};
return $points->[$#$points] if @$points;
}
## @method NumPoints()
# @brief Return the number of points in the sequence.
#
# @note Returns all points in a list context.
sub NumPoints {
my($self) = @_;
@{$self->{Points}};
}
## @method PointN($N, $point)
# @param N A location in the sequence
# @note The first point has the index 1 as OGC SF SQL conformance test uses 1-based indexing.
# @param point [optional] A Point object, if defined sets the point to index N
sub PointN {
my($self, $N, $point) = @_;
my $points = $self->{Points};
$points->[$N-1] = $point if defined $point;
return $points->[$N-1];
}
sub Is3D {
my($self) = @_;
for my $p (@{$self->{Points}}) {
return 1 if $p->Is3D;
}
return 0;
}
sub IsMeasured {
my($self) = @_;
for my $p (@{$self->{Points}}) {
return 1 if $p->IsMeasured;
}
return 0;
}
sub IsClosed {
my($self) = @_;
$self->StartPoint()->Equals($self->EndPoint());
}
## @method Close()
# @brief Close the curve by adding the first point also as the last point.
# @note Not in the specification.
sub Close {
my($self) = @_;
push @{$self->{Points}}, $self->{Points}[0];
}
## @method IsRing($upgrade)
# @brief Tests whether this curve is a ring, i.e., closed and simple
# @param upgrade [optional, not in the specification] Upgrades this
# curve into a Ring if this really could be a ring
sub IsRing {
my($self, $upgrade) = @_;
my $ret = ($self->IsClosed and $self->IsSimple);
bless($self, 'Geo::OGC::LinearRing') if $ret and $upgrade;
return $ret;
}
# should use Precision if one exists
sub Equals {
my($self, $geom) = @_;
return 0 unless $geom->isa('Geo::OGC::Curve');
return 0 unless $#{$self->{Points}} == $#{$geom->{Points}};
for my $i (0..$#{$self->{Points}}) {
return 0 unless $self->{Points}[$i]->Equals($geom->{Points}[$i]);
}
return 1;
}
sub Area {
return 0;
}
sub MakeCollection {
my($self) = @_;
my $coll = Geo::OGC::MultiCurve->new;
$coll->AddGeometry($self);
return $coll;
}
sub ApplyTransformation {
my($self, $transf) = @_;
for my $p (@{$self->{Points}}) {
$p->ApplyTransformation($transf);
}
}
## @method Reverse()
# @brief Reverse the order of the points in the sequence.
# @note Not in the specification.
sub Reverse {
my($self) = @_;
@{$self->{Points}} = reverse @{$self->{Points}};
}
sub ClosestVertex {
my($self, $x, $y) = @_;
return unless @{$self->{Points}};
my($dmin) = $self->{Points}[0]->ClosestVertex($x, $y);
my $i = 0;
for my $j (1..$#{$self->{Points}}) {
my($d) = $self->{Points}[$j]->ClosestVertex($x, $y);
($i, $dmin) = ($j, $d) if $d < $dmin;
}
return ($i, $dmin);
}
sub VertexAt {
my($self, $i) = @_;
return ($self->{Points}[0], $self->{Points}[$#{$self->{Points}}])
if (($i == 0 or $i == $#{$self->{Points}}) and $self->isa('Geo::OGC::LinearRing'));
return ($self->{Points}[$i]);
}
sub _closest_point {
my($x0, $y0, $x1, $y1, $x, $y) = @_;
my $ab2 = ($x1-$x0)*($x1-$x0) + ($y1-$y0)*($y1-$y0);
my $ap_ab = ($x-$x0)*($x1-$x0) + ($y-$y0)*($y1-$y0);
my $t = $ap_ab/$ab2;
if ($t < 0) {$t = 0} elsif ($t > 1) {$t = 1}
my $xp = $x0+$t*($x1-$x0);
my $yp = $y0+$t*($y1-$y0);
return ($xp, $yp, ($x-$xp)*($x-$xp)+($y-$yp)*($y-$yp));
}
sub ClosestPoint {
my($self, $x, $y) = @_;
return unless @{$self->{Points}};
my($i, $pmin, $dmin);
for my $j (1..$#{$self->{Points}}) {
my($xp, $yp, $d) =
_closest_point($self->{Points}[$j-1]{X}, $self->{Points}[$j-1]{Y},
$self->{Points}[$j]{X}, $self->{Points}[$j]{Y}, $x, $y);
($i, $pmin, $dmin) = ($j, Geo::OGC::Point->new($xp, $yp), $d)
if (!defined($dmin) or $d < $dmin);
}
return ($i, $pmin, $dmin)
}
sub AddVertex {
my($self, $i, $p) = @_;
splice @{$self->{Points}}, $i, 0, $p;
}
sub DeleteVertex {
my($self, $i) = @_;
splice @{$self->{Points}}, $i, 1;
}
#
# LineString
#
package Geo::OGC::LineString;
use strict;
use Carp;
use Geo::OGC::Geometry qw/:all/;
our @ISA = qw( Geo::OGC::Curve );
# de Berg et al p. 25
sub FindIntersections {
# @_ contains a list of linestrings that make up the S
my @linestrings = @_;
#my $precision =
# event queue
my @Q = (); # [ymax,index1,index2], ..., index1 is index to @_ index2 is index to Points
for my $index1 (0 .. $#linestrings) {
my $s = $linestrings[$index1]->{Points};
for my $index2 (0 .. $#$s-1) {
my($y1, $y2) = ( $s->[$index2]{Y}, $s->[$index2+1]{Y} );
if ($y1 > $y2) {
push @Q, [$y1, $index1, $index2];
push @Q, [$y2];
} else {
push @Q, [$y2, $index1, $index2];
push @Q, [$y1];
}
}
}
# process event points in descending ymax order
@Q = sort {$b->[0] <=> $a->[0]} @Q;
my $T = Tree::Binary::Search->new(); # the status structure
$T->setComparisonFunction
(sub {
my($a, $b) = @_; # two keys
});
my $i = 0;
while ($i < @Q) {
my $j = $i+1;
while ($j < @Q and sqrt(($Q[$i][0]-$Q[$j][0])**2) < $SNAP_DISTANCE_SQR) {
$j++;
}
# $i .. $j-1 are the event points in @Q
for my $k ($i .. $j-1) {
if (@{$Q[$k]} == 3) {
my $i1 = $Q[$k][1];
my $i2 = $Q[$k][2];
my $s = $linestrings[$i1]->{Points};
$T->insert( "$i1,$i2" => [$s->[$i2]{X}, $s->[$i2]{Y}, $s->[$i2+1]{X}, $s->[$i2+1]{Y}] );
}
}
}
}
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::Curve::new($package, %params);
return $self;
}
sub GeometryType {
return 'LineString';
}
sub IsSimple {
my($self) = @_;
my $edges = @{$self->{Points}} - 1;
return 1 if $edges < 2;
my $closed = $self->IsClosed;
my $simple = 1;
for my $i (0..$edges-1) {
# check a case of self tangency
return 0 if $i < $edges-1 and $self->{Points}[$i+2]->Equals($self->{Points}[$i]);
for my $j ($i+2..$edges-1) {
next if $closed and $i == 0 and $j == $edges-1;
return 0 if intersect
(
$self->{Points}[$i]{X}, $self->{Points}[$i]{Y},
$self->{Points}[$i+1]{X}, $self->{Points}[$i+1]{Y},
$self->{Points}[$j]{X}, $self->{Points}[$j]{Y},
$self->{Points}[$j+1]{X},$self->{Points}[$j+1]{Y}
);
}
}
return 1;
}
sub Envelope {
my($self) = @_;
my($minx, $miny, $maxx, $maxy);
for my $p (@{$self->{Points}}) {
$minx = $p->{X} if !defined($minx) or $minx > $p->{X};
$miny = $p->{Y} if !defined($miny) or $miny > $p->{Y};
$maxx = $p->{X} if !defined($maxx) or $maxx > $p->{X};
$maxy = $p->{Y} if !defined($maxy) or $maxy > $p->{Y};
}
my $r = Geo::OGC::LinearRing->new;
$r->AddPoint(Geo::OGC::Point->new($minx, $miny));
$r->AddPoint(Geo::OGC::Point->new($maxx, $miny));
$r->AddPoint(Geo::OGC::Point->new($maxx, $maxy));
$r->AddPoint(Geo::OGC::Point->new($minx, $maxy));
$r->Close;
return $r;
}
## @method Length()
# @brief The length of this LineString in its associated spatial reference.
# @note Currently computed as a simple euclidean distance.
sub Length {
my($self) = @_;
my $l = 0;
my($x0, $y0) = ($self->{Points}[0]{X}, $self->{Points}[0]{Y});
for my $i (1..$#{$self->{Points}}) {
my($x1, $y1) = ($self->{Points}[$i]{X}, $self->{Points}[$i]{Y});
$l += sqrt(($x1 - $x0)**2+($y1 - $y0)**2);
($x0, $y0) = ($x1, $y1);
}
return $l;
}
sub Distance {
my($self, $geom) = @_;
if ($geom->isa('Geo::OGC::Point')) {
return $geom->DistanceToLineStringSqr($self);
} elsif ($geom->isa('Geo::OGC::LineString')) {
my $dist;
for my $p (@{$self->{Points}}) {
my $d = $p->DistanceToLineStringSqr($geom);
$dist = $d if !(defined $dist) or $d < $dist;
}
return $dist;
} elsif ($geom->isa('Geo::OGC::Polygon')) {
my $dist;
for my $p (@{$self->{Points}}) {
my $d = $p->Distance($geom);
$dist = $d if !(defined $dist) or $d < $dist;
}
return $dist;
} elsif ($geom->isa('Geo::OGC::GeometryCollection')) {
my $dist = $self->Distance($geom->{Geometries}[0]);
for my $g (@{$geom->{Geometries}}[1..$#{$geom->{Geometries}}]) {
my $d = $self->Distance($g);
$dist = $d if $d < $dist;
}
return $dist;
} else {
croak "can't compute distance between a ".ref($geom)." and a line string";
}
}
sub LinesWhereWithin {
my($self, $point) = @_;
my($x, $y) = ($point->{X}, $point->{Y});
my @ret;
my $p1 = $self->{Points}[0];
return @ret unless $p1;
my $p2 = $self->{Points}[1];
return @ret unless $p1;
push @ret, 1 if
distance_point_line_sqr($x, $y, $p1->{X}, $p1->{Y}, $p2->{X}, $p2->{Y}) < $SNAP_DISTANCE_SQR;
for my $i (2..$#{$self->{Points}}) {
$p1 = $p2;
$p2 = $self->{Points}[$i];
push @ret, $i if
distance_point_line_sqr($x, $y, $p1->{X}, $p1->{Y}, $p2->{X}, $p2->{Y}) < $SNAP_DISTANCE_SQR;
}
return @ret;
}
sub Within {
my($self, $geom) = @_;
if ($geom->isa('Geo::OGC::Point')) {
for my $p (@{$self->{Points}}) {
return 0 unless $p->Equals($geom);
}
return 1;
} elsif ($geom->isa('Geo::OGC::LineString')) {
my @w1 = ();
for my $p (@{$self->{Points}}) {
my @w2 = $geom->LinesWhereWithin($p);
return 0 unless @w2;
next unless @w1;
my $overlap = 0;
for my $w1 (@w1) {
for my $w2 (@w2) {
$overlap = 1 if $w1 == $w2;
}
}
return 0 unless $overlap;
@w1 = @w2;
}
return 1;
} elsif ($geom->isa('Geo::OGC::Polygon')) {
for my $p (@{$self->{Points}}) {
return 0 if $geom->{ExteriorRing}->IsPointStricktlyOut($p);
for my $ring (@{$geom->{InteriorRings}}) {
return 0 if $ring->IsPointStricktlyIn($p);
}
}
my $i = $self->Intersection($geom->{ExteriorRing});
for my $g (@{$i->{Geometries}}) {
next unless $g->isa('Geo::OGC::Line');
# does the line go out of the polygon?
# yes, if its start and end points are on different lines
my @s = $geom->{ExteriorRing}->LinesWhereWithin($g->StartPoint);
my @e = $geom->{ExteriorRing}->LinesWhereWithin($g->EndPoint);
my $overlap = 0;
for my $s (@s) {
for my $e (@e) {
$overlap = 1 if $s == $e;
}
}
return 0 unless $overlap;
}
for my $ring (@{$geom->{InteriorRings}}) {
my $i = $self->Intersection($ring);
for my $g (@{$i->{Geometries}}) {
next unless $g->isa('Geo::OGC::Line');
# does the line go into the interior ring?
# yes, if its start and end points are on different lines
my @s = $ring->LinesWhereWithin($g->StartPoint);
my @e = $ring->LinesWhereWithin($g->EndPoint);
my $overlap = 0;
for my $s (@s) {
for my $e (@e) {
$overlap = 1 if $s == $e;
}
}
return 0 unless $overlap;
}
}
return 1;
} else {
croak "linestring within ".ref($geom)." is not yet implemented";
}
}
# assuming simple lines
# z and m coords!
sub Intersection {
my($self, $geom) = @_;
if ($geom->isa('Geo::OGC::Point')) {
return $geom->Clone if $geom->DistanceToLineStringSqr($self) < $SNAP_DISTANCE_SQR;
} elsif ($geom->isa('Geo::OGC::LineString')) {
#my $i = Geo::OGC::GeometryCollection->new;
my %i;
my $index = 1;
my $p1;
for my $p2 (@{$self->{Points}}) {
$p1 = $p2, next unless $p1;
my $q1;
for my $q2 (@{$geom->{Points}}) {
$q1 = $q2, next unless $q1;
my @p = ($p1->{X}, $p1->{Y}, $p2->{X}, $p2->{Y});
my @q = ($q1->{X}, $q1->{Y}, $q2->{X}, $q2->{Y});
#print STDERR "lines: @p; @q\n";
if (intersect( @p, @q )) {
my $p1q = distance_point_line_sqr(@p[0..1], @q);
my $p2q = distance_point_line_sqr(@p[2..3], @q);
my $q1p = distance_point_line_sqr(@q[0..1], @p);
my $q2p = distance_point_line_sqr(@q[2..3], @p);
#print STDERR "p1=q $p1q p2=q $p2q q1=p $q1p q2=p $q2p\n";
if ($p1q < $SNAP_DISTANCE_SQR) {
if ($q1p < $SNAP_DISTANCE_SQR) {
if ($p1->Equals($q1)) {
$i{$index++} = $p1->Clone;
} else {
$i{$index++} = Geo::OGC::Line->new( points=>[[@p[0..1]],[@q[0..1]]] );
}
}
if ($q2p < $SNAP_DISTANCE_SQR) {
if ($p1->Equals($q2)) {
$i{$index++} = $p1->Clone;
} else {
$i{$index++} = Geo::OGC::Line->new( points=>[[@p[0..1]],[@q[2..3]]] );
}
}
if ($p2q < $SNAP_DISTANCE_SQR) {
$i{$index++} = Geo::OGC::Line->new( points=>[[@p[0..1]],[@p[2..3]]] );
} else {
$i{$index++} = Geo::OGC::Point->new( @p[0..1] );
}
} elsif ($p2q < $SNAP_DISTANCE_SQR) {
if ($q1p < $SNAP_DISTANCE_SQR) {
if ($p2->Equals($q1)) {
$i{$index++} = $p2->Clone;
} else {
$i{$index++} = Geo::OGC::Line->new( points=>[[@p[2..3]],[@q[0..1]]] );
}
}
if ($q2p < $SNAP_DISTANCE_SQR) {
if ($p2->Equals($q2)) {
$i{$index++} = $p2->Clone;
} else {
$i{$index++} = Geo::OGC::Line->new( points=>[[@p[2..3]],[@q[2..3]]] );
}
}
$i{$index++} = Geo::OGC::Point->new( @p[2..3] );
} elsif ($q1p < $SNAP_DISTANCE_SQR) {
if ($q2p < $SNAP_DISTANCE_SQR) {
$i{$index++} = Geo::OGC::Line->new( points=>[[@q[0..1]],[@q[2..3]]] );
} else {
$i{$index++} = Geo::OGC::Point->new( @q[0..1] );
}
} elsif ($q2p < $SNAP_DISTANCE_SQR) {
$i{$index++} = Geo::OGC::Point->new( @q[2..3] );
} else {
$i{$index++} = Geo::OGC::Point->new( Geo::OGC::Geometry::intersection_point(@p, @q) );
}
}
$q1 = $q2;
}
$p1 = $p2;
}
#$i->Simplify;
# delete unnecessary points and lines
# comparisons are done unnecessarily twice??
for my $g1 (keys %i) {
for my $g2 (keys %i) {
next if $g1 == $g2;
if ($i{$g1}->Within($i{$g2})) {
#print STDERR "delete ",$i{$g1}->AsText," because it is in ",$i{$g2}->AsText,"\n";
delete $i{$g1};
last;
}
}
}
my $i = Geo::OGC::GeometryCollection->new;
for my $g1 (keys %i) {
$i->AddGeometry($i{$g1});
}
return $i;
} else {
croak "intersection between a ".ref($geom)." and a line string is not yet implemented";
}
}
sub MakeCollection {
my($self) = @_;
my $coll = Geo::OGC::MultiLineString->new;
$coll->AddGeometry($self);
return $coll;
}
# from http://everything2.com/index.pl?node_id=859282
sub pt_to_seg_dist {
# distance of p to segment p1-p2, v12 is the vector p1p2
my ($p1, $v12, $p) = @_;
my $m12 = $v12->[0] * $v12->[0] + $v12->[1] * $v12->[1];
my $v1p = [];
$v1p->[0] = $p->{X} - $p1->{X};
$v1p->[1] = $p->{Y} - $p1->{Y};
my $dot = $v1p->[0] * $v12->[0] + $v1p->[1] * $v12->[1];
if ($dot <= 0.0)
{
return sqrt ($v1p->[0] * $v1p->[0] + $v1p->[1] * $v1p->[1]);
}
else
{
if ($dot >= $m12)
{
$v1p->[0] = $v1p->[0] + $v12->[0];
$v1p->[1] = $v1p->[1] + $v12->[1];
return sqrt ($v1p->[0] * $v1p->[0] + $v1p->[1] * $v1p->[1]);
}
else
{
my $slash = $v1p->[0] * $v12->[1] - $v1p->[1] * $v12->[0];
return abs ($slash / sqrt ($m12));
}
}
}
sub simplify_part {
my($self, $first, $last, $simple, $tolerance) = @_;
if ($last > $first + 1)
{
my $p1 = $self->{Points}[$first];
my $vfl = [$self->{Points}[$last]{X} - $self->{Points}[$first]{X},
$self->{Points}[$last]{Y} - $self->{Points}[$first]{Y}];
# find the intermediate point
# furthest from the segment
# connecting first and last
my $b = $first+1;
my $db = pt_to_seg_dist ($p1, $vfl, $self->{Points}[$b]);
my $i = $b + 1;
while ($i < $last)
{
my $di = pt_to_seg_dist ($p1, $vfl, $self->{Points}[$i]);
if ($di > $db)
{
$b = $i;
$db = $di;
}
$i++;
}
# if the furthest distance beats the tolerance,
# recursively simplify the rest of the array.
if ($db >= $tolerance)
{
simplify_part ($self, $first, $b, $simple, $tolerance);
$simple->AddPoint($self->{Points}[$b]);
simplify_part ($self, $b, $last, $simple, $tolerance);
}
}
}
## @method simplify($tolerance)
# @brief Simplifies the linestring using Douglas-Peucker
# @return The simpliefied linestring
sub simplify {
my($self, $tolerance) = @_;
my $simple = Geo::OGC::LineString->new;
$simple->AddPoint($self->StartPoint);
simplify_part ($self, 0, $self->NumPoints-1, $simple, $tolerance);
$simple->AddPoint($self->EndPoint);
return $simple;
}
#
# Line
#
package Geo::OGC::Line;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::LineString );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::LineString::new($package, %params);
return $self;
}
sub GeometryType {
return 'Line';
}
#
# LinearRing
#
package Geo::OGC::LinearRing;
use strict;
use Carp;
use Geo::OGC::Geometry qw/:all/;
our @ISA = qw( Geo::OGC::LineString );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::LineString::new($package, %params);
return $self;
}
sub GeometryType {
return 'LinearRing';
}
## @method IsPointIn(Geo::OGC::Point point)
# @brief Tests whether the given point is within the ring
# @note uses the pnpoly algorithm from
# http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
# @note Assumes a simple closed ring
# @note may or may not return true if the point is on the border
sub IsPointIn {
my($self, $point) = @_;
my($x, $y) = ($point->{X}, $point->{Y});
my $c = 0;
my $prev;
for my $p (@{$self->{Points}}) {
$prev = $p, next unless $prev;
$c = !$c if (((( $p->{Y} <= $y ) && ( $y < $prev->{Y} )) ||
(( $prev->{Y} <= $y ) && ( $y < $p->{Y} ))) &&
( $x < ( $prev->{X} - $p->{X} ) *
( $y - $p->{Y} ) / ( $prev->{Y} - $p->{Y} ) + $p->{X} ));
$prev = $p;
}
return $c;
}
# @note not on the border
sub IsPointStricktlyIn {
my($self, $point) = @_;
return 1 if ( $self->IsPointIn($point) and
!($point->DistanceToLineStringSqr($self) < $SNAP_DISTANCE_SQR) );
return 0;
}
# @note not on the border
sub IsPointStricktlyOut {
my($self, $point) = @_;
return 1 unless ( $self->IsPointIn($point) or
$point->DistanceToLineStringSqr($self) < $SNAP_DISTANCE_SQR );
return 0;
}
## @method Area()
# @brief Computes the area of the ring.
# @note Not in the specification.
# @note Computed as a simple euclidean area.
# @note Assumes a simple closed ring
# @return area or negative area if the sense of the rotation of the ring is clockwise
sub Area {
my($self) = @_;
my $area = 0;
my $N = $#{$self->{Points}}-1; # skip the closing point
my $j = 0;
for my $i (0..$N) {
$j++;
$area += $self->{Points}[$i]{X} * $self->{Points}[$j]{Y};
$area -= $self->{Points}[$i]{Y} * $self->{Points}[$j]{X};
}
return $area/2;
}
sub Centroid {
my($self) = @_;
my($area, $x, $y) = (0, 0, 0);
my $N = $#{$self->{Points}}-1; # skip the closing point
for my $i (0..$N) {
my($xi, $yi, $xj, $yj) = ( $self->{Points}[$i]{X}, $self->{Points}[$i]{Y},
$self->{Points}[$i+1]{X}, $self->{Points}[$i+1]{Y} );
my $b = $xi * $yj - $xj * $yi;
$area += $b;
$x += ($xi + $xj) * $b;
$y += ($yi + $yj) * $b;
}
$x /= abs(3*$area); # 6 but $area is 2 * real area
$y /= abs(3*$area);
return Geo::OGC::Point->new($x, $y);
}
## @method IsCCW()
# @brief Returns true if the points in this ring are arranged counterclockwise.
# @note Not in the specification.
# @note Assumes a simple closed ring.
sub IsCCW {
my($self) = @_;
# find the northernmost point
my $t = 0;
my $N = $#{$self->{Points}}-1; # skip the closing point
for my $i (1..$N) {
$t = $i if $self->{Points}[$i]{Y} > $self->{Points}[$t]{Y};
}
# the previous point
my $p = $t-1;
$p = $N if $p < 0;
# the next point
my $n = $t+1;
$n = 0 if $n > $N;
return ccw($self->{Points}[$p]{X}, $self->{Points}[$p]{Y},
$self->{Points}[$t]{X}, $self->{Points}[$t]{Y},
$self->{Points}[$n]{X}, $self->{Points}[$n]{Y}) == 1;
}
## @method Rotate()
# @brief Makes clockwise from counterclockwise and vice versa.
sub Rotate {
my($self) = @_;
@{$self->{Points}} = reverse @{$self->{Points}};
}
#
# Surface
#
package Geo::OGC::Surface;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::Geometry );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::Geometry::new($package, %params);
return $self;
}
sub GeometryType {
return 'Surface';
}
sub Dimension {
return 2;
}
sub Area {
my($self) = @_;
croak "Area method for class ".ref($self)." is not implemented yet";
}
sub Centroid {
my($self) = @_;
croak "Centroid method for class ".ref($self)." is not implemented yet";
}
sub PointOnSurface {
my($self) = @_;
croak "PointOnSurface method for class ".ref($self)." is not implemented yet";
}
sub MakeCollection {
my($self) = @_;
my $coll = Geo::OGC::MultiSurface->new;
$coll->AddGeometry($self);
return $coll;
}
#
# Polygon
#
package Geo::OGC::Polygon;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::Surface );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::Surface::new($package, %params);
return $self;
}
sub init {
my($self, %params) = @_;
$self->SUPER::init(%params);
$self->{ExteriorRing} = undef;
$self->{InteriorRings} = [];
}
sub copy {
my($self, $clone) = @_;
$self->SUPER::copy($clone);
$clone->ExteriorRing($self->{ExteriorRing}->Clone) if $self->{ExteriorRing};
for my $r (@{$self->{InteriorRings}}) {
$clone->AddInteriorRing($r->Clone);
}
}
sub GeometryType {
return 'Polygon';
}
## @method Assert()
# @brief Test the rules that define valid polygons
# @note Not in the specification.
sub Assert {
my($self) = @_;
# a) Polygons are topologically closed;
croak "not at least triangle" unless @{$self->{ExteriorRing}->{Points}} > 3;
croak "exterior not closed" unless $self->{ExteriorRing}->IsClosed;
# b) The boundary of a Polygon consists of a set of LinearRings
# that make up its exterior and interior boundaries;
for my $r (@{$self->{InteriorRings}}) {
croak "an interior is not closed" unless $r->IsClosed;
}
# c) No two Rings in the boundary cross and the Rings in the
# boundary of a Polygon may intersect at a Point but only as a
# tangent
for my $ring (@{$self->{InteriorRings}}) {
for my $p (@{$ring->{Points}}) {
croak "point in interior not within exterior" unless $self->{ExteriorRing}->IsPointIn($p);
for my $r2 (@{$self->{InteriorRings}}) {
next if $ring == $r2;
croak "point in interior is within another interior" if $r2->IsPointIn($p);
}
}
}
# d) A Polygon may not have cut lines, spikes or punctures
croak "exterior is not simple" unless $self->{ExteriorRing}->IsSimple;
for my $r (@{$self->{InteriorRings}}) {
croak "an interior is not simple" unless $r->IsSimple;
}
# e) The interior of every Polygon is a connected point set
# f) The exterior of a Polygon with 1 or more holes is not
# connected. Each hole defines a connected component of the
# exterior.
for my $i (0..$#{$self->{InteriorRings}}) {
my $r1 = $self->{InteriorRings}[$i];
my $r2 = $r1->Intersection($self->{ExteriorRing});
croak "an interior intersects too much with the exterior"
if @{$r2->{Geometries}} and (@{$r2->{Geometries}} > 1 or
!$r2->{Geometries}[0]->isa('Geo::OGC::Point'));
for my $j ($i+1..$#{$self->{InteriorRings}}) {
my $r2 = $self->{InteriorRings}[$j];
$r2 = $r1->Intersection($r2);
croak "an interior intersects too much with another interior"
if @{$r2->{Geometries}} and (@{$r2->{Geometries}} > 1 or
!$r2->{Geometries}[0]->isa('Geo::OGC::Point'));
}
}
return 1;
}
sub Is3D {
my($self) = @_;
return 1 if $self->{ExteriorRing}->Is3D;
for my $r (@{$self->{InteriorRings}}) {
return 1 if $r->Is3D;
}
return 0;
}
sub IsMeasured {
my($self) = @_;
return 1 if $self->{ExteriorRing}->IsMeasured;
for my $r (@{$self->{InteriorRings}}) {
return 1 if $r->IsMeasured;
}
return 0;
}
sub AddInteriorRing {
my($self, $ring, $i) = @_;
croak 'usage: Polygon->AddInteriorRing($ring[, $i])'
unless $ring and $ring->isa('Geo::OGC::LinearRing');
my $rings = $self->{InteriorRings};
$i = @$rings unless defined $i;
if (@$rings) {
my $temp = $rings->[$i-1];
splice @$rings,$i-1,1,($temp, $ring);
} else {
push @$rings, $ring;
}
}
sub ExteriorRing {
my($self, $ring) = @_;
if (defined $ring) {
croak 'usage: Polygon->ExteriorRing($ring)'
unless $ring->isa('Geo::OGC::LinearRing');
$self->{ExteriorRing} = $ring;
} else {
return $self->{ExteriorRing};
}
}
sub Envelope {
my($self) = @_;
return $self->{ExteriorRing}->Envelope;
}
## @method NumInteriorRing()
# @brief Return the number of interior rings in the polygon.
#
# @note Returns all interior rings in a list context.
sub NumInteriorRing {
my($self) = @_;
@{$self->{InteriorRings}};
}
sub InteriorRingN {
my($self, $N, $ring) = @_;
my $rings = $self->{InteriorRings};
$rings->[$N-1] = $ring if defined $ring;
return $rings->[$N-1] if @$rings;
}
# @note Assumes the order of the interior rings is the same.
sub Equals {
my($self, $geom) = @_;
return 0 unless $geom->isa('Geo::OGC::Polygon');
return 0 unless @{$self->{InteriorRings}} == @{$geom->{InteriorRings}};
return 0 unless $self->{ExteriorRing}->Equals($geom->{ExteriorRing});
for my $i (0..$#{$self->{InteriorRings}}) {
return 0 unless $self->{InteriorRings}[$i]->Equals($geom->{InteriorRings}[$i]);
}
return 1;
}
sub Area {
my($self) = @_;
my $a = $self->{ExteriorRing}->Area;
for my $ring (@{$self->{InteriorRings}}) {
$a -= $ring->Area;
}
return $a;
}
sub IsPointIn {
my($self, $point) = @_;
my $c = $self->{ExteriorRing}->IsPointIn($point);
if ($c) {
for my $ring (@{$self->{InteriorRings}}) {
$c = 0, last if $ring->IsPointIn($point);
}
}
return $c;
}
sub Distance {
my($self, $geom) = @_;
if ($geom->isa('Geo::OGC::Point')) {
return $geom->Distance($self);
} elsif ($geom->isa('Geo::OGC::LineString')) {
return $geom->Distance($self);
} elsif ($geom->isa('Geo::OGC::Polygon')) {
my $dist;
for my $p (@{$self->{ExteriorRing}->{Points}}) {
if ($geom->{ExteriorRing}->IsPointIn($p)) {
my $c = 1;
for my $ring (@{$self->{InteriorRings}}) {
if ($ring->IsPointIn($p)) {
my $d = $p->DistanceToLineStringSqr($ring);
$dist = $d if !(defined $dist) or $d < $dist;
$c = 0;
}
}
return 0 if $c;
} else {
my $d = $p->DistanceToLineStringSqr($geom->{ExteriorRing});
$dist = $d if !(defined $dist) or $d < $dist;
}
}
return $dist;
} elsif ($geom->isa('Geo::OGC::GeometryCollection')) {
my $dist = $self->Distance($geom->{Geometries}[0]);
for my $g (@{$geom->{Geometries}}[1..$#{$geom->{Geometries}}]) {
my $d = $self->Distance($g);
$dist = $d if $d < $dist;
}
return $dist;
} else {
croak "can't compute distance between a ".ref($geom)." and a polygon";
}
}
sub Within {
my($self, $geom) = @_;
if ($geom->isa('Geo::OGC::Point')) {
for my $p (@{$self->{ExteriorRing}->{Points}}) {
return 0 unless $p->Equals($geom);
}
return 1;
} elsif ($geom->isa('Geo::OGC::LineString')) {
for my $p (@{$self->{ExteriorRing}->{Points}}) {
return 0 unless $p->Within($geom);
}
return 1;
} elsif ($geom->isa('Geo::OGC::Polygon')) {
croak "polygon within ".ref($geom)." is not yet implemented";
# the exterior and interior rings must be completely within
# and the other's interiors must not be within ...
for my $p (@{$self->{ExteriorRing}->{Points}}) {
return 0 unless $p->Within($geom);
}
return 1;
} else {
croak "polygon within ".ref($geom)." is not yet implemented";
}
}
sub Intersection {
my($self, $geom) = @_;
if ($geom->isa('Geo::OGC::Point')) {
return $geom->Intersection($self);
} elsif ($geom->isa('Geo::OGC::LineString')) {
return $geom->Intersection($self);
} elsif ($geom->isa('Geo::OGC::Polygon')) {
croak "polygon within ".ref($geom)." is not yet implemented";
# 1st intersection between outer rings $A (this) and $B (other)
my $A = $self->{ExteriorRing};
my $B = $geom->{ExteriorRing};
my $p1 = $A->{Points}[0];
my $c = $B->IsPointStricktlyIn($p1);
my $i = 0;
my $j;
while (1) {
my $p = $A->{Points}[$i];
$j = $i, last if $B->IsPointStricktlyOut($p);
$i++;
last if $i == @{$A->{Points}};
}
if (defined $j) { # there is at least one point
}
# the list of
my @A; # lines on A
} else {
croak "polygon within ".ref($geom)." is not yet implemented";
}
}
sub as_text {
my($self, $force_parens, $include_tag) = @_;
my $text .= $self->{ExteriorRing}->as_text if $self->{ExteriorRing};
for my $ring (@{$self->{InteriorRings}}) {
$text .= ', ';
$text .= $ring->as_text;
}
my $ZM = $self->Is3D ? 'Z' : '';
$ZM .= 'M' if $self->IsMeasured;
$text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
return $text;
}
sub MakeCollection {
my($self) = @_;
my $coll = Geo::OGC::MultiPolygon->new;
$coll->AddGeometry($self);
return $coll;
}
sub ClosestVertex {
my($self, $x, $y) = @_;
return unless $self->{ExteriorRing};
my($imin, $dmin) = $self->{ExteriorRing}->ClosestVertex($x, $y);
my $iring = -1;
my $r = 0;
for my $ring (@{$self->{InteriorRings}}) {
my($i, $d) = $ring->ClosestVertex($x, $y);
($iring, $imin, $dmin) = ($r, $i, $d) if $d < $dmin;
$r++;
}
return ($iring, $imin, $dmin);
}
sub VertexAt {
my($self, $iring, $ivertex) = @_;
return $self->{ExteriorRing}->VertexAt($ivertex) if $iring < 0;
return $self->{InteriorRings}[$iring]->VertexAt($ivertex);
}
sub ClosestPoint {
my($self, $x, $y) = @_;
return unless $self->{ExteriorRing};
my($imin, $pmin, $dmin) = $self->{ExteriorRing}->ClosestPoint($x, $y);
my $iring = -1;
my $r = 0;
for my $ring (@{$self->{InteriorRings}}) {
my($i, $p, $d) = $ring->ClosestPoint($x, $y);
($iring, $imin, $pmin, $dmin) = ($r, $i, $p, $d) if $d < $dmin;
$r++;
}
return ($iring, $imin, $pmin, $dmin);
}
sub AddVertex {
my($self, $ring, $i, $p) = @_;
$self->{ExteriorRing}->AddVertex($i, $p), return if $ring < 0;
$self->{InteriorRings}[$ring]->AddVertex($i, $p);
}
sub DeleteVertex {
my($self, $ring, $i) = @_;
$self->{ExteriorRing}->DeleteVertex($i), return if $ring < 0;
$self->{InteriorRings}[$ring]->DeleteVertex($i);
}
## @method LastPolygon()
# @brief Returns self
sub LastPolygon {
my($self) = @_;
return $self;
}
#
# Triangle
#
package Geo::OGC::Triangle;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::Polygon );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::Polygon::new($package, %params);
return $self;
}
sub GeometryType {
return 'Triangle';
}
#
# PolyhedralSurface
#
package Geo::OGC::PolyhedralSurface;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::Surface );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::Surface::new($package, %params);
return $self;
}
sub init {
my($self, %params) = @_;
$self->SUPER::init(%params);
$self->{Patches} = []; # polygon
if ($params{patches}) {
for my $p (@{$params{patches}}) {
$self->AddPatch(Geo::OGC::Polygon->new(points=>$p));
}
} elsif ($params{patchesm}) {
for my $p (@{$params{patches}}) {
$self->AddPatch(Geo::OGC::Polygon->new(pointsm=>$p));
}
}
}
sub copy {
my($self, $clone) = @_;
$self->SUPER::copy($clone);
for my $p (@{$self->{Patches}}){
$clone->AddPatch($p->Clone);
}
}
sub GeometryType {
return 'PolyhedralSurface';
}
sub AddPatch {
my($self, $patch, $i) = @_;
croak 'usage: PolyhedralSurface->AddPatch($patch[, $i])'
unless $patch and $patch->isa('Geo::OGC::Polygon');
my $patches = $self->{Patches};
$i = @$patches unless defined $i;
if (@$patches) {
my $temp = $patches->[$i-1];
splice @$patches,$i-1,1,($temp, $patch);
} else {
push @$patches, $patch;
}
}
sub NumPatches {
my($self) = @_;
@{$self->{Patches}};
}
sub PatchN {
my($self, $N, $patch) = @_;
my $patches = $self->{Patches};
$patches->[$N-1] = $patch if defined $patch;
return $patches->[$N-1] if @$patches;
}
sub BoundingPolygons {
my($self, $p) = @_;
croak "BoundingPolygons method for class ".ref($self)." is not implemented yet";
}
sub IsClosed {
my($self) = @_;
croak "IsClosed method for class ".ref($self)." is not implemented yet";
}
sub IsMeasured {
my($self) = @_;
for my $p (@{$self->{Patches}}) {
return 1 if $p->IsMeasured;
}
return 0;
}
sub as_text {
my($self, $force_parens, $include_tag) = @_;
my $text = '';
for my $patch (@{$self->{Patches}}) {
$text .= ', ';
$text .= $patch->as_text;
}
my $ZM = $self->Is3D ? 'Z' : '';
$ZM .= 'M' if $self->IsMeasured;
$text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
return $text;
}
#
# TIN
#
package Geo::OGC::TIN;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::PolyhedralSurface );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::Surface::new($package, %params);
return $self;
}
sub GeometryType {
return 'TIN';
}
#
# GeometryCollection
#
package Geo::OGC::GeometryCollection;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::Geometry );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::Geometry::new($package, %params);
return $self;
}
sub init {
my($self, %params) = @_;
$self->SUPER::init(%params);
$self->{Geometries} = [];
}
sub copy {
my($self, $clone) = @_;
$self->SUPER::copy($clone);
for my $g (@{$self->{Geometries}}) {
$clone->AddGeometry($g->Clone);
}
}
sub GeometryType {
return 'GeometryCollection';
}
sub Dimension {
my($self) = @_;
my $dim;
for my $g (@{$self->{Geometries}}) {
my $d = $g->Dimension;
$dim = $d if !(defined $dim) or $d > $dim;
}
return $dim;
}
sub Is3D {
my($self) = @_;
for my $g (@{$self->{Geometries}}) {
return 1 if $g->Is3D;
}
return 0;
}
sub IsMeasured {
my($self) = @_;
for my $g (@{$self->{Geometries}}) {
return 1 if $g->IsMeasured;
}
return 0;
}
sub as_text {
my($self, $force_parens, $include_tag) = @_;
my $text = join(',', map {$_->as_text(1, 1)} @{$self->{Geometries}});
my $ZM = $self->Is3D ? 'Z' : '';
$ZM .= 'M' if $self->IsMeasured;
$text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
return $text;
}
sub ElementType {
return 'Geometry';
}
sub AddGeometry {
my($self, $geometry, $i) = @_;
croak 'usage: GeometryCollection->AddGeometry($geometry[, $i])'
unless $geometry and $geometry->isa('Geo::OGC::Geometry');
my $geometries = $self->{Geometries};
$i = @$geometries unless defined $i;
if (@$geometries) {
my $temp = $geometries->[$i-1];
splice @$geometries,$i-1,1,($temp, $geometry);
} else {
push @$geometries, $geometry;
}
}
## @method NumGeometries()
# @brief Return the number of geometries in the collection.
#
# @note Returns all geometries in a list context.
sub NumGeometries {
my($self) = @_;
@{$self->{Geometries}};
}
## @method GeometryN(N)
# @brief Return the Nth geometry from the collection
# (the index of the first geometry is 1).
#
sub GeometryN {
my($self, $N, $geometry) = @_;
my $geometries = $self->{Geometries};
$geometries->[$N-1] = $geometry if defined $geometry;
return $geometries->[$N-1] if @$geometries;
}
sub Envelope {
my($self) = @_;
my($minx, $miny, $maxx, $maxy);
for my $g (@{$self->{Geometries}}) {
my $e = $g->Envelope;
my $min = $e->PointN(1);
my $max = $e->PointN(3);
$minx = $min->{X} if !defined($minx) or $minx > $min->{X};
$miny = $min->{Y} if !defined($miny) or $miny > $min->{Y};
$maxx = $max->{X} if !defined($maxx) or $maxx > $max->{X};
$maxy = $max->{Y} if !defined($maxy) or $maxy > $max->{Y};
}
my $r = new Geo::OGC::LinearRing;
$r->AddPoint(Geo::OGC::Point->new($minx, $miny));
$r->AddPoint(Geo::OGC::Point->new($maxx, $miny));
$r->AddPoint(Geo::OGC::Point->new($maxx, $maxy));
$r->AddPoint(Geo::OGC::Point->new($minx, $maxy));
$r->Close;
return $r;
}
# @note Assumes the order is the same.
sub Equals {
my($self, $geom) = @_;
return 0 unless $geom->isa('Geo::OGC::GeometryCollection');
return 0 unless @{$self->{Geometries}} == @{$geom->{Geometries}};
for my $i (0..$#{$self->{Geometries}}) {
return 0 unless $self->{Geometries}[$i]->Equals($geom->{Geometries}[$i]);
}
return 1;
}
sub Distance {
my($self, $geom) = @_;
if ($geom->isa('Geo::OGC::Point')) {
return $geom->Distance($self);
} elsif ($geom->isa('Geo::OGC::LineString')) {
return $geom->Distance($self);
} elsif ($geom->isa('Geo::OGC::Polygon')) {
return $geom->Distance($self);
} elsif ($geom->isa('Geo::OGC::GeometryCollection')) {
my $dist = $self->Distance($geom->{Geometries}[0]);
for my $g (@{$geom->{Geometries}}[1..$#{$geom->{Geometries}}]) {
my $d = $g->Distance($self);
$dist = $d if $d < $dist;
}
return $dist;
} else {
croak "can't compute distance between a ".ref($geom)." and a geometry collection";
}
}
sub ClosestVertex {
my($self, $x, $y) = @_;
return unless @{$self->{Geometries}};
my @rmin = $self->{Geometries}[0]->ClosestVertex($x, $y);
my $imin = 0;
my $i = 0;
for my $g (@{$self->{Geometries}}) {
$i++, next if $i == 0;
my @r = $g->ClosestVertex($x, $y);
if ($r[$#r] < $rmin[$#rmin]) {
@rmin = @r;
$imin = $i;
}
$i++;
}
return ($imin, @rmin);
}
sub VertexAt {
my $self = shift;
my $i = shift;
return $self->{Geometries}[$i]->VertexAt(@_);
}
sub ClosestPoint {
my($self, $x, $y) = @_;
return unless @{$self->{Geometries}};
my @rmin = $self->{Geometries}[0]->ClosestPoint($x, $y);
my $imin = 0;
my $i = 0;
for my $g (@{$self->{Geometries}}) {
$i++, next if $i == 0;
my @r = $g->ClosestPoint($x, $y);
if ($r[$#r] < $rmin[$#rmin]) {
@rmin = @r;
$imin = $i;
}
$i++;
}
return ($imin, @rmin);
}
sub AddVertex {
my $self = shift;
my $i = shift;
$self->{Geometries}[$i]->AddVertex(@_);
}
sub DeleteVertex {
my $self = shift;
my $i = shift;
$self->{Geometries}[$i]->DeleteVertex(@_);
}
## @method LastPolygon()
# @brief Returns last polygon or undef
sub LastPolygon {
my($self) = @_;
for (my $i = $#{$self->{Geometries}}; $i >= 0; $i--) {
my $polygon = $self->{Geometries}[$i]->LastPolygon;
return $polygon if $polygon;
}
}
#
# MultiSurface
#
package Geo::OGC::MultiSurface;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::GeometryCollection );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::GeometryCollection::new($package, %params);
return $self;
}
sub GeometryType {
return 'MultiSurface';
}
sub ElementType {
return 'Surface';
}
sub Area {
my($self) = @_;
croak "Area method for class ".ref($self)." is not implemented yet";
}
sub Centroid {
my($self) = @_;
croak "Centroid method for class ".ref($self)." is not implemented yet";
}
sub PointOnSurface {
my($self) = @_;
croak "PointOnSurface method for class ".ref($self)." is not implemented yet";
}
#
# MultiCurve
#
package Geo::OGC::MultiCurve;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::GeometryCollection );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::GeometryCollection::new($package, %params);
return $self;
}
sub GeometryType {
return 'MultiCurve';
}
sub ElementType {
return 'Curve';
}
#
# MultiPoint
#
package Geo::OGC::MultiPoint;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::GeometryCollection );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::GeometryCollection::new($package, %params);
return $self;
}
sub init {
my($self, %params) = @_;
$self->SUPER::init(%params);
if ($params{points}) {
for my $p (@{$params{points}}) {
$self->AddGeometry(Geo::OGC::Point->new(point=>$p));
}
} elsif ($params{pointsm}) {
for my $p (@{$params{pointsm}}) {
$self->AddGeometry(Geo::OGC::Point->new(pointm=>$p));
}
}
}
sub GeometryType {
return 'MultiPoint';
}
sub as_text {
my($self, $force_parens, $include_tag) = @_;
my $text = join(',', map {$_->as_text(1)} @{$self->{Geometries}});
my $ZM = $self->Is3D ? 'Z' : '';
$ZM .= 'M' if $self->IsMeasured;
$text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
return $text;
}
sub ElementType {
return 'Point';
}
sub Boundary {
my($self) = @_;
return Geo::OGC::GeometryCollection->new();
}
#
# MultiPolygon
#
package Geo::OGC::MultiPolygon;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::MultiSurface );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::MultiSurface::new($package, %params);
return $self;
}
sub GeometryType {
return 'MultiPolygon';
}
sub as_text {
my($self, $force_parens, $include_tag) = @_;
my $text = join(',', map {$_->as_text(1)} @{$self->{Geometries}});
my $ZM = $self->Is3D ? 'Z' : '';
$ZM .= 'M' if $self->IsMeasured;
$text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
return $text;
}
sub ElementType {
return 'Polygon';
}
#
# MultiLineString
#
package Geo::OGC::MultiLineString;
use strict;
use Carp;
our @ISA = qw( Geo::OGC::MultiCurve );
sub new {
my($package, %params) = @_;
my $self = Geo::OGC::MultiCurve::new($package, %params);
return $self;
}
sub GeometryType {
return 'MultiLineString';
}
sub as_text {
my($self, $force_parens, $include_tag) = @_;
my $text = join(',', map {$_->as_text(1)} @{$self->{Geometries}});
my $ZM = $self->Is3D ? 'Z' : '';
$ZM .= 'M' if $self->IsMeasured;
$text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
return $text;
}
sub ElementType {
return 'LineString';
}