Skip to content

Commit

Permalink
Merge pull request #13 from OpenSmock/dev
Browse files Browse the repository at this point in the history
  • Loading branch information
@labordep
labordep authored Feb 29, 2024
2 parents 9d5e320 + bf4efdc commit 8adec9f
Show file tree
Hide file tree
Showing 4 changed files with 133 additions and 10 deletions.
2 changes: 1 addition & 1 deletion GeoView-Bloc-Alexandrie/GeoViewUserToDisplayToGraphicLayer.extension.st
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,5 +7,5 @@ GeoViewUserToDisplayToGraphicLayer >> aeDrawOn: aeCanvas [
sortedDatas ifNil:[ ^ self ].

self flag: 'Patch to wait correction on issue #9'.
sortedDatas copy do: [ :e | e aeDrawOn: aeCanvas ]
sortedDatas copy do: [ :e | e ifNotNil:[:el | el aeDrawOn: aeCanvas ]]
]
87 changes: 87 additions & 0 deletions GeoView-Tests/GeodesicUtilitiesTest.class.st
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,87 @@
"
A GeodesicUtilitiesTest is a test class for testing the behavior of GeodesicUtilities
"
Class {
#name : #GeodesicUtilitiesTest,
#superclass : #TestCase,
#category : #'GeoView-Tests-MapProjection'
}

{ #category : #tests }
GeodesicUtilitiesTest >> testConvertGeodesicToAzimuthInRadiansFromTwoPoints [

| originPointAbsCoord endPointAbsCoord azRad|

originPointAbsCoord := AbsoluteCoordinates latitudeInDegrees: 48.28275 longitudeInDegrees: -4.70553.
endPointAbsCoord := AbsoluteCoordinates latitudeInDegrees: 48.32676 longitudeInDegrees: -4.68284.

azRad := GeodesicUtilities convertGeodesicToAzimuthInRadiansFrom: originPointAbsCoord to: endPointAbsCoord.

self assert: (azRad radiansToDegrees closeTo: 18.9726921431 precision: 1e-10) equals: true.

endPointAbsCoord := AbsoluteCoordinates latitudeInDegrees: 48.29603 longitudeInDegrees: -4.62295.

azRad := GeodesicUtilities convertGeodesicToAzimuthInRadiansFrom: originPointAbsCoord to: endPointAbsCoord.

self assert: (azRad radiansToDegrees closeTo: 76.42078534394 precision: 1e-10) equals: true.
]

{ #category : #tests }
GeodesicUtilitiesTest >> testConvertGeodesicToDistanceFromTwoPoints [

| originPointAbsCoord endPointAbsCoord distanceInM|

originPointAbsCoord := AbsoluteCoordinates latitudeInDegrees: 48.28275 longitudeInDegrees: -4.70553.
endPointAbsCoord := AbsoluteCoordinates latitudeInDegrees: 48.32676 longitudeInDegrees: -4.68284.

distanceInM := GeodesicUtilities convertGeodesicToDistanceInMeterFrom: originPointAbsCoord to: endPointAbsCoord.

self assert: (distanceInM closeTo: 5175.1422410892 precision: 1e-8) equals: true.

endPointAbsCoord := AbsoluteCoordinates latitudeInDegrees: 48.29603 longitudeInDegrees: -4.62295.

distanceInM := GeodesicUtilities convertGeodesicToDistanceInMeterFrom: originPointAbsCoord to: endPointAbsCoord.

self assert: (distanceInM closeTo: 6303.42663738019 precision: 1e-8) equals: true.
]

{ #category : #tests }
GeodesicUtilitiesTest >> testIsAzimuthTowardsArea [

| originPointAbsCoord isAzTowardsArea listPointsArea headingInDegree|

originPointAbsCoord := AbsoluteCoordinates latitudeInDegrees: 48.28275 longitudeInDegrees: -4.70553.

listPointsArea := OrderedCollection new.
listPointsArea add: (AbsoluteCoordinates latitudeInDegrees: 48.32676 longitudeInDegrees: -4.68284).
listPointsArea add: (AbsoluteCoordinates latitudeInDegrees: 48.33995 longitudeInDegrees: -4.65976).
listPointsArea add: (AbsoluteCoordinates latitudeInDegrees: 48.3117 longitudeInDegrees: -4.60466).
listPointsArea add: (AbsoluteCoordinates latitudeInDegrees: 48.29603 longitudeInDegrees: -4.62295).

"headingInDegree := 45.
isAzTowardsArea := GeodesicUtilities isAzimuthTowardsAreaFrom: originPointAbsCoord azimuth: headingInDegree area: listPointsArea.
self assert: isAzTowardsArea equals: true.
originPointAbsCoord latitudeInDegrees: 48.29603.
isAzTowardsArea := GeodesicUtilities isAzimuthTowardsAreaFrom: originPointAbsCoord azimuth: headingInDegree area: listPointsArea.
self assert: isAzTowardsArea equals: true.
headingInDegree := -45.
isAzTowardsArea := GeodesicUtilities isAzimuthTowardsAreaFrom: originPointAbsCoord azimuth: headingInDegree area: listPointsArea.
self assert: isAzTowardsArea equals: false."
originPointAbsCoord latitudeInDegrees: 48.29603.
headingInDegree := 90.

isAzTowardsArea := GeodesicUtilities isAzimuthTowardsAreaFrom: originPointAbsCoord azimuth: headingInDegree area: listPointsArea.

self assert: isAzTowardsArea equals: true.



]
48 changes: 39 additions & 9 deletions GeoView/GeodesicUtilities.class.st
View file
Original file line number Diff line number Diff line change
Expand Up @@ -153,12 +153,12 @@ GeodesicUtilities class >> convertGeodesicToAzimuthInRadiansFrom: startAbsoluteC
(startAbsoluteCoordinate longitudeInRadians closeTo: endAbsoluteCoordinate longitudeInRadians precision: 1e-7) ifTrue:
[ (startAbsoluteCoordinate latitudeInRadians closeTo: endAbsoluteCoordinate latitudeInRadians precision: 1e-7)
ifTrue: [^ 0].
endAbsoluteCoordinate latitudeInRadians >= startAbsoluteCoordinate latitudeInRadians abs
endAbsoluteCoordinate latitudeInRadians >= startAbsoluteCoordinate latitudeInRadians
ifTrue: [^ 0]
ifFalse: [^ Float pi]].

d5 := ((startAbsoluteCoordinate latitudeInRadians tan * WGS84 semiMinorAxisInMeters) / WGS84 semiMajorAxisInMeters) arcTan.
d6 := ((startAbsoluteCoordinate longitudeInRadians tan * WGS84 semiMinorAxisInMeters) / WGS84 semiMajorAxisInMeters) arcTan.
d6 := ((endAbsoluteCoordinate latitudeInRadians tan * WGS84 semiMinorAxisInMeters) / WGS84 semiMajorAxisInMeters) arcTan.
d7 := d5 cos.
d8 := d6 cos.
d9 := d5 sin.
Expand All @@ -177,7 +177,8 @@ GeodesicUtilities class >> convertGeodesicToAzimuthInRadiansFrom: startAbsoluteC
d17 := (d7 * d10) - (d9 * d8 * d15).
d20 := ((d16 * d16) + (d17 * d17)) sqrt.
d21 := (d9 * d10) + (d7 * d8 * d15).
d22 := (d20 @ d21) theta.
" !!! Math.atan2(y,x) => (x,y) theta "
d22 := (d21 @ d20) theta.
d23 := (d7 * d8 * d14) / d20.
d24 := 1 - (d23 * d23).
d25 := d24 < 1e-14
Expand All @@ -191,9 +192,10 @@ GeodesicUtilities class >> convertGeodesicToAzimuthInRadiansFrom: startAbsoluteC
d14 := d12 sin.
d15 := d12 cos.
break := (d12 - d19) abs < 1e-10.
d28 := d28 + 1 ].
break ifFalse:[d28 := d28 + 1 ]
].

d28 := ((d8 * d14) @ ((d7 * d10) - (d9 * d8 * d15))) theta.
d28 := ( ((d7 * d10) - (d9 * d8 * d15)) @ (d8 * d14)) theta.
^ d28 negative
ifTrue: [ d28 + (2 * Float pi) ]
ifFalse: [ d28 ]
Expand All @@ -205,7 +207,11 @@ GeodesicUtilities class >> convertGeodesicToDistanceInMeterFrom: startAbsoluteCo
"Compute distance between two Geodesic (lat/long) coordinates and return in meter"
<tacplot60>
| d5 d6 d7 d8 d9 d10 d11 d12 d14 d15 d16 d17 d18 d19 d20 d21 d22 d23 d24 d25 d26 d27 d28 d29 d30 d31 d32 xxx l |
startAbsoluteCoordinate = endAbsoluteCoordinate ifTrue: [ ^0 ].

(startAbsoluteCoordinate longitudeInRadians closeTo: endAbsoluteCoordinate longitudeInRadians precision: 1e-7) ifTrue:
[ (startAbsoluteCoordinate latitudeInRadians closeTo: endAbsoluteCoordinate latitudeInRadians precision: 1e-7)
ifTrue: [^ 0]
].

d5 := ((startAbsoluteCoordinate latitudeInRadians tan * WGS84 semiMinorAxisInMeters) / WGS84 semiMajorAxisInMeters) arcTan.
d6 := ((endAbsoluteCoordinate latitudeInRadians tan * WGS84 semiMinorAxisInMeters) / WGS84 semiMajorAxisInMeters) arcTan.
Expand All @@ -222,13 +228,14 @@ GeodesicUtilities class >> convertGeodesicToDistanceInMeterFrom: startAbsoluteCo

l := 1.
xxx := 1.
[(xxx > 1e-10 and: [ l < 4 ])] whileTrue:
[(xxx >= 1e-10 and: [ l < 4 ])] whileTrue:
[ xxx := (d12 - d20) abs.
d16 := d8 * d14.
d18 := (d7 * d10) - (d9 * d8 * d15).
d21 := ((d16 * d16) + (d18 * d18)) sqrt.
d22 := (d9 * d10) + (d7 * d8 * d15).
d23 := (d21 @ d22) theta.
" !!! Math.atan2(y,x) => (x,y) theta "
d23 := (d22 @ d21) theta.
d24 := (d7 * d8 * d14) / d21.
d25 := 1 - (d24 * d24).
d26 := d25 < 1e-14
Expand All @@ -243,7 +250,7 @@ GeodesicUtilities class >> convertGeodesicToDistanceInMeterFrom: startAbsoluteCo
d15 := d12 cos.
l := l + 1. ].

d29 := d25 * (WGS84 semiMajorAxisInMeters - WGS84 semiMinorAxisInMeters) / WGS84 semiMinorAxisInMeters.
d29 := d25 * WGS84 eminorSquare.
d30 := 1 + ((d29 / 16384) * (4096 + (d29 * (-768 + (d29 * (320 - (175 * d29))))))).
d31 := (d29 / 1024) * (256 + (d29 * (-128 + (d29 * (74 - (47 * d29)))))).
d17 := d22 * (-1 + (2 * (d26 * d26))).
Expand All @@ -252,6 +259,29 @@ GeodesicUtilities class >> convertGeodesicToDistanceInMeterFrom: startAbsoluteCo
^ WGS84 semiMinorAxisInMeters * d30 * (d23 - d32)
]

{ #category : #tools }
GeodesicUtilities class >> isAzimuthTowardsAreaFrom: anAbsCoordOrigin azimuth: anAzimuthInDegree area: aListOfAbsoluteCoords [

| listAngles isPointingArea azRad |
isPointingArea := false.

aListOfAbsoluteCoords isEmptyOrNil ifTrue: [ ^ isPointingArea ].

"Compute angles between origin position and vertices of area"
listAngles := (aListOfAbsoluteCoords collect: [ :absCoord |
self convertGeodesicToAzimuthInRadiansFrom: anAbsCoordOrigin to: absCoord ]
) asOrderedCollection.

listAngles := listAngles sort: [ :a :b | a < b ].

azRad := anAzimuthInDegree degreesToRadians.
isPointingArea := azRad between: listAngles first and: listAngles last.
isPointingArea := isPointingArea or:[ azRad closeTo: listAngles first precision: 1e-3 ].
isPointingArea := isPointingArea or:[ azRad closeTo: listAngles last precision: 1e-3 ].

^ isPointingArea
]

{ #category : #tools }
GeodesicUtilities class >> localPolar: aPolarCoordinates toGeodesicFrom: anOriginAbsoluteCoordinates [
"Compute the geodesic position of a specified line end, knowing its azimuth, elevation, length and start geodesic position. Warning: that method doesn't take in account the earth curve as orthodromic2geodesic does.
Expand Down
6 changes: 6 additions & 0 deletions GeoView/WGS84.class.st
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,12 @@ Class {
#category : #'GeoView-MapProjection'
}

{ #category : #'standard values' }
WGS84 class >> eminorSquare [

^ ((self semiMajorAxisInMeters * self semiMajorAxisInMeters) - (self semiMinorAxisInMeters * self semiMinorAxisInMeters)) / ((self semiMinorAxisInMeters * self semiMinorAxisInMeters)).
]

{ #category : #'standard values' }
WGS84 class >> flattening [
"https://en.wikipedia.org/wiki/World_Geodetic_System"
Expand Down

0 comments on commit 8adec9f

Please sign in to comment.