source: subversion/applications/utils/gary68/mwCoastLines.pm @ 34714

Last change on this file since 34714 was 26259, checked in by gary68, 8 years ago

new mapweaver version

File size: 11.2 KB
Line 
1#
2# PERL mapweaver module by gary68
3#
4#
5#
6#
7# Copyright (C) 2011, Gerhard Schwanz
8#
9# 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
10# Free Software Foundation; either version 3 of the License, or (at your option) any later version.
11#
12# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License along with this program; if not, see <http://www.gnu.org/licenses/>
16#
17
18
19package mwCoastLines ; 
20
21use strict ;
22use warnings ;
23use Math::Polygon ;
24use List::Util qw[min max] ;
25
26use mwMap ;
27use mwFile ;
28use mwConfig ;
29use mwMisc ;
30
31use vars qw($VERSION @ISA @EXPORT @EXPORT_OK);
32
33require Exporter ;
34
35@ISA = qw ( Exporter AutoLoader ) ;
36
37@EXPORT = qw (  processCoastLines
38
39                 ) ;
40
41
42sub nearestPoint {
43#
44# accepts x/y coordinates and returns nearest point on border of map to complete cut coast ways
45#
46        my $ref = shift ;
47        my $x = $ref->[0] ;
48        my $y = $ref->[1] ;
49        my $xn ; my $yn ;
50        my $min = 99999 ;
51        # print "  NP: initial $x $y\n" ;
52        my ($xmax, $ymax) = getDimensions() ;
53        # print "  NP: dimensions $xmax $ymax\n" ;
54        if ( abs ($xmax-$x) < $min) { # right
55                $xn = $xmax ;
56                $yn = $y ; 
57                $min = abs ($xmax-$x) ;
58        }
59        if ( abs ($ymax-$y) < $min) { # bottom
60                $xn = $x ;
61                $yn = $ymax ; 
62                $min = abs ($ymax-$y) ;
63        }
64        if ( abs ($x) < $min) { # left
65                $xn = 0 ;
66                $yn = $y ; 
67                $min = abs ($x) ;
68        }
69        if ( abs ($y) < $min) { # top
70                $xn = $x ;
71                $yn = 0 ; 
72        }
73        # print "  NP: final $xn $yn\n" ;
74        my @a = ($xn, $yn) ;
75        return (\@a) ;
76}
77
78
79
80sub nextPointOnBorder {
81#
82# accepts x/y coordinates and returns next point on border - to complete coast rings with other polygons and corner points
83# hints if returned point is a corner
84#
85        # right turns
86        my ($x, $y) = @_ ;
87        my ($xn, $yn) ;
88        my $corner = 0 ;
89        my ($xmax, $ymax) = getDimensions() ;
90        if ($x == $xmax) { # right border
91                if ($y < $ymax) {
92                        $xn = $xmax ; $yn = $y + 1 ;
93                }
94                else {
95                        $xn = $xmax - 1 ; $yn = $ymax ;
96                }
97        }
98        else {
99                if ($x == 0) { # left border
100                        if ($y > 0) {
101                                $xn = 0 ; $yn = $y - 1 ;
102                        }
103                        else {
104                                $xn = 1 ; $yn = 0 ;
105                        }
106                }
107                else {
108                        if ($y == $ymax) { # bottom border
109                                if ($x > 0) {
110                                        $xn = $x - 1 ; $yn = $ymax ;
111                                }
112                                else {
113                                        $xn = 0 ; $yn = $ymax - 1 ; 
114                                }
115                        }
116                        else {
117                                if ($y == 0) { # top border
118                                        if ($x < $xmax) {
119                                                $xn = $x + 1 ; $yn = 0 ;
120                                        }
121                                        else {
122                                                $xn = $xmax ; $yn = 1 ; 
123                                        }
124                                }
125                        }
126                }
127        }
128        # print "NPOB: $x, $y --- finito $xn $yn\n" ;
129
130        if ( ($xn == 0) and ($yn == 0) ) { $corner = 1 ; }
131        if ( ($xn == 0) and ($yn == $ymax) ) { $corner = 1 ; }
132        if ( ($xn == $xmax) and ($yn == 0) ) { $corner = 1 ; }
133        if ( ($xn == $xmax) and ($yn == $ymax) ) { $corner = 1 ; }
134
135        return ($xn, $yn, $corner) ;
136}
137
138# ---------------------------------------------------------------------------------
139
140sub processCoastLines {
141#
142#
143#
144        print "check and process coastlines...\n" ;
145
146        my $ref = shift ; # ref to all coast ways
147        my @allWays = @$ref ;
148
149        if (cv('debug') eq "1") { 
150                print "COAST: " . scalar (@allWays) . " coast ways initially found.\n" ;
151                print "COAST: ways: @allWays\n\n" ;
152        }
153
154
155        my ($lonRef, $latRef) = getNodePointers() ;
156        my ($nodesRef, $tagRef) = getWayPointers() ;
157
158        # check coast ways. eliminate invisible ways. eliminate points outside map.
159        my @newWays = () ;
160        foreach my $w ( @allWays ) {
161                my @nodes = @{ $$nodesRef{ $w } } ;
162
163                my $allIn = 1 ;
164                my $allOut = 1 ; 
165                foreach my $n ( @nodes ) {
166                        if ( pointInMap ($n) ) {
167                                $allOut = 0 ;
168                        }
169                        else {
170                                $allIn = 0 ;
171                        }
172                }
173
174                if ( $allIn ) {
175                        # use way as it is     
176                        push @newWays, $w ;
177                        if ( cv ('debug') eq "1" ) { print "COAST: way $w will be used unmodified.\n" ; }
178                }
179                elsif ( $allOut) {
180                        # do nothing
181                        if ( cv ('debug') eq "1" ) { print "COAST: way $w will NOT be used. outside map.\n" ; }
182                }
183                else {
184                        # eliminate all outside nodes at start and end of way, then use new way
185                       
186                        # eliminate outsides at start
187                        while ( (scalar @nodes >= 1) and ( ! pointInMap ($nodes[0]) ) ) {
188                                shift @nodes ;
189                        }
190
191                        # eliminate outsides at end
192                        while ( (scalar @nodes >= 1) and ( ! pointInMap ($nodes[-1]) ) ) {
193                                pop @nodes ;
194                        }
195
196                        if ( scalar @nodes >= 2 ) {
197                                @{ $$nodesRef{$w}} = @nodes ;
198                                push @newWays, $w ;
199                                if ( cv ('debug') eq "1" ) { print "COAST: modified way $w will be used.\n" ; }
200                        }
201                        else {
202                                if ( cv ('debug') eq "1" ) { print "COAST: way $w too short now.\n" ; }
203                        }
204
205                }               
206
207        }
208
209        @allWays = @newWays ;
210
211
212
213        if (cv('debug') eq "1") { 
214                print "\nCOAST: " . scalar (@allWays) . " coast ways will be used.\n" ;
215                print "COAST: ways: @allWays\n\n" ;
216        }
217
218        if (scalar @allWays > 0) {
219                # build rings
220                my ($refWays, $refNodes) = buildRings (\@allWays, 0) ;
221                my @ringNodes = @$refNodes ; # contains all nodes of rings // array of arrays !
222                if (cv('debug') eq "1") { print "COAST: " . scalar (@ringNodes) . " rings found.\n" ; }
223
224                # convert rings to coordinate system
225                my @ringCoordsOpen = () ; my @ringCoordsClosed = () ;
226                for (my $i=0; $i<=$#ringNodes; $i++) {
227                        # print "COAST: initial ring $i\n" ;
228                        my @actualCoords = () ;
229                        foreach my $node (@{$ringNodes[$i]}) {
230                                push @actualCoords, [convert ($$lonRef{$node}, $$latRef{$node})] ;
231                        }
232                        if (${$ringNodes[$i]}[0] == ${$ringNodes[$i]}[-1]) {
233                                push @ringCoordsClosed, [@actualCoords] ; # islands
234                        }
235                        else {
236                                push @ringCoordsOpen, [@actualCoords] ;
237                        }
238                        # printRingCoords (\@actualCoords) ;
239                        my $num = scalar @actualCoords ;
240                        if (cv('debug') eq "1") { print "COAST: initial ring $i - $actualCoords[0]->[0],$actualCoords[0]->[1] -->> $actualCoords[-1]->[0],$actualCoords[-1]->[1]  nodes: $num\n" ; }
241                }
242
243                if (cv('debug') eq "1") { print "COAST: add points on border...\n" ; }
244                foreach my $ring (@ringCoordsOpen) {
245                        # print "COAST:   ring $ring with border nodes\n" ;
246                        # add first point on border
247                        my $ref = nearestPoint ($ring->[0]) ;
248                        my @a = @$ref ;
249                        unshift @$ring, [@a] ;
250                        # add last point on border
251                        $ref = nearestPoint ($ring->[-1]) ;
252                        @a = @$ref ;
253                        push @$ring, [@a] ;
254                        # printRingCoords ($ring) ;
255                }
256
257                my @islandRings = @ringCoordsClosed ;
258                if (cv('debug') eq "1") { print "COAST: " . scalar (@islandRings) . " islands found.\n" ; }
259                @ringCoordsClosed = () ;
260
261                # process ringCoordsOpen
262                # add other rings, corners...
263                while (scalar @ringCoordsOpen > 0) { # as long as there are open rings
264                        if (cv('debug') eq "1") { print "COAST: building ring...\n" ; }
265                        my $ref = shift @ringCoordsOpen ; # get start ring
266                        my @actualRing = @$ref ;
267
268                        my $closed = 0 ; # mark as not closed
269                        my $actualX = $actualRing[-1]->[0] ;
270                        my $actualY = $actualRing[-1]->[1] ;
271
272                        my $actualStartX = $actualRing[0]->[0] ; 
273                        my $actualStartY = $actualRing[0]->[1] ; 
274
275                        if (cv('debug') eq "1") { print "COAST: actual and actualStart $actualX, $actualY   -   $actualStartX, $actualStartY\n" ; }
276
277                        my $corner ;
278                        while (!$closed) { # as long as this ring is not closed
279                                ($actualX, $actualY, $corner) = nextPointOnBorder ($actualX, $actualY) ;
280                                # print "      actual $actualX, $actualY\n" ;
281                                my $startFromOtherPolygon = -1 ;
282                                # find matching ring if there is another ring
283                                if (scalar @ringCoordsOpen > 0) {
284                                        for (my $i=0; $i <= $#ringCoordsOpen; $i++) {
285                                                my @test = @{$ringCoordsOpen[$i]} ;
286                                                # print "    test ring $i: ", $test[0]->[0], " " , $test[0]->[1] , "\n" ;
287                                                if ( ($actualX == $test[0]->[0]) and ($actualY == $test[0]->[1]) ) {
288                                                        $startFromOtherPolygon = $i ;
289                                                        if (cv('debug') eq "1") { print "COAST:   matching start other polygon found i= $i\n" ; }
290                                                }
291                                        }
292                                }
293                                # process matching polygon, if present
294                                if ($startFromOtherPolygon != -1) { # start from other polygon {
295                                        # append nodes
296                                        # print "ARRAY TO PUSH: @{$ringCoordsOpen[$startFromOtherPolygon]}\n" ;
297                                        push @actualRing, @{$ringCoordsOpen[$startFromOtherPolygon]} ;
298                                        # set actual
299                                        $actualX = $actualRing[-1]->[0] ;
300                                        $actualY = $actualRing[-1]->[1] ;
301                                        # drop p2 from opens
302                                        splice @ringCoordsOpen, $startFromOtherPolygon, 1 ;
303                                        if (cv('debug') eq "1") { print "COAST:   openring $startFromOtherPolygon added to actual ring\n" ; }
304                                }
305                                else {
306                                        if ($corner) { # add corner to actual ring
307                                                push @actualRing, [$actualX, $actualY] ;
308                                                if (cv('debug') eq "1") { print "COAST:   corner $actualX, $actualY added to actual ring\n" ; }
309                                        }
310                                }
311                                # check if closed
312                                if ( ($actualX == $actualStartX) and ($actualY == $actualStartY) ) {
313                                        $closed = 1 ;
314                                        push @actualRing, [$actualX, $actualY] ;
315                                        push @ringCoordsClosed, [@actualRing] ;
316                                        if (cv('debug') eq "1") { print "COAST:    ring now closed and moved to closed rings.\n" ; }
317                                }
318                        } # !closed
319                } # open rings
320
321                my $color = cv('oceancolor') ;
322
323                # build islandRings polygons
324                if (cv('debug') eq "1") { print "OCEAN: building island polygons\n" ; }
325                my @islandPolygons = () ;
326                if (scalar @islandRings > 0) {
327                        for (my $i=0; $i<=$#islandRings; $i++) {
328                                my @poly = () ;
329                                foreach my $node ( @{$islandRings[$i]} ) {
330                                        push @poly, [$node->[0], $node->[1]] ;
331                                }
332                                my ($p) = Math::Polygon->new(@poly) ;
333                                $islandPolygons[$i] = $p ;
334                        }
335                }
336               
337                # build ocean ring polygons
338                if (cv('debug') eq "1") { print "OCEAN: building ocean polygons\n" ; }
339                my @oceanPolygons = () ;
340                if (scalar @ringCoordsClosed > 0) {
341                        for (my $i=0; $i<=$#ringCoordsClosed; $i++) {
342                                my @poly = () ;
343                                foreach my $node ( @{$ringCoordsClosed[$i]} ) {
344                                        push @poly, [$node->[0], $node->[1]] ;
345                                }
346                                my ($p) = Math::Polygon->new(@poly) ;
347                                $oceanPolygons[$i] = $p ;
348                        }
349                }
350                else {
351                        if (scalar @islandRings > 0) {
352                                if (cv('debug') eq "1") { print "OCEAN: build ocean rect\n" ; }
353                                my @ocean = () ;
354                                my ($x, $y) = getDimensions() ;
355                                push @ocean, [0,0], [$x,0], [$x,$y], [0,$y], [0,0] ;
356                                push @ringCoordsClosed, [@ocean] ;
357                                my ($p) = Math::Polygon->new(@ocean) ;
358                                push @oceanPolygons, $p ;
359                        }
360                }
361
362                # finally create pathes for SVG
363                for (my $i=0; $i<=$#ringCoordsClosed; $i++) {
364                # foreach my $ring (@ringCoordsClosed) {
365                        my @ring = @{$ringCoordsClosed[$i]} ;
366                        my @array = () ;
367                        my @coords = () ;
368                        foreach my $c (@ring) {
369                                push @coords, $c->[0], $c->[1] ;
370                        }
371                        push @array, [@coords] ; 
372                        if (scalar @islandRings > 0) {
373                                for (my $j=0; $j<=$#islandRings; $j++) {
374                                        # island in ring? 1:1 and coast on border?
375                                        # if (isIn ($islandPolygons[$j], $oceanPolygons[$i]) == 1) {
376                                        if ( (isIn ($islandPolygons[$j], $oceanPolygons[$i]) == 1) or 
377                                                ( (scalar @islandRings == 1) and (scalar @ringCoordsClosed == 1) ) )    {
378                                                if (cv('debug') eq "1") { print "OCEAN: island $j in ocean $i\n" ; }
379                                                my @coords = () ;
380                                                foreach my $c (@{$islandRings[$j]}) {
381                                                        push @coords, $c->[0], $c->[1] ;               
382                                                }
383                                                push @array, [@coords] ;
384                                        }
385                                }
386                        }
387
388
389                        # drawAreaOcean ($color, \@array) ;
390                        my $svgText = "fill=\"$color\" " ;
391                        drawArea($svgText, "none", \@array, 0, "base") ;
392
393                }
394        }
395}
396
397sub pointInMap {
398        my ($n) = shift ;
399        my ($sizeX, $sizeY) = getDimensions() ;
400        my ($lonRef, $latRef) = getNodePointers() ;
401
402        my ($x, $y) = convert ($$lonRef{$n}, $$latRef{$n}) ;
403
404        my $ok = 0 ;
405        if (
406                ( $x >= 0 ) and
407                ( $x <= $sizeX ) and
408                ( $y >= 0 ) and
409                ( $y <= $sizeY ) ) {
410                $ok = 1 ;
411        }
412        return $ok ;           
413}
414
4151 ;
416
417
Note: See TracBrowser for help on using the repository browser.