Because of numerical floating point issues, we round the points to a …

…specified number of significant digits when doing checks
GeoRegionsEcosystem · Feb 14, 2025 · 6aaab4b · 6aaab4b
1 parent 482da97
commit 6aaab4b
Show file tree

Hide file tree

Showing 2 changed files with 54 additions and 22 deletions.
diff --git a/src/isin/isin.jl b/src/isin/isin.jl
@@ -2,7 +2,8 @@
     in(
         point :: Point2{<:Real},
         geo   :: GeoRegion;
-        throw :: Bool = false
+        throw :: Bool = false,
+        sigdigits :: Int = 10
     ) -> tf :: Bool
 
 Check if a geographical point `point` is within a GeoRegion defined by `geo`.
@@ -15,6 +16,7 @@ Arguments
 Keyword Arguments
 =================
 - `throw` : If `true`, then if `point` is not within `geo`, an error is thrown and the program stops running.
+- `sigdigits` : Specifies number of significant digits (i.e., precision) of the point coordinates used for checking. Defaults to 10.
 
 Returns
 =======
@@ -23,24 +25,28 @@ Returns
 function Base.in(
     point :: Point2{<:Real},
     geo   :: GeoRegion;
-    throw :: Bool = false
+    throw :: Bool = false,
+    sigdigits :: Int = 10
 )
 
     throw ? (@info "$(modulelog()) - Performing a check to determine if the coordinates $(point) are within the specified region boundaries.") : nothing
 
-    plon = point[1]
-    plat = point[2]
+    plon = point[1]; plat = round(point[2],sigdigits=sigdigits)
 
     while plon > 360;  plon -= 360 end
     while plon < -180; plon += 360 end
 
+    plon1 = round(plon,sigdigits=sigdigits)
+    plon2 = round(plon-360,sigdigits=sigdigits)
+    plon3 = round(plon+360,sigdigits=sigdigits)
+
     isin = !iszero(sum([
-        within(Point(plon    ,plat),geo.geometry.polygon),
-        within(Point(plon+360,plat),geo.geometry.polygon),
-        within(Point(plon-360,plat),geo.geometry.polygon),
-        touches(Point(plon    ,plat),geo.geometry.polygon),
-        touches(Point(plon+360,plat),geo.geometry.polygon),
-        touches(Point(plon-360,plat),geo.geometry.polygon)
+        within(Point(plon1,plat),geo.geometry.polygon),
+        within(Point(plon2,plat),geo.geometry.polygon),
+        within(Point(plon3,plat),geo.geometry.polygon),
+        touches(Point(plon1,plat),geo.geometry.polygon),
+        touches(Point(plon2,plat),geo.geometry.polygon),
+        touches(Point(plon3,plat),geo.geometry.polygon)
     ]))
 
     if !isin
@@ -65,7 +71,8 @@ end
         geo  :: GeoRegion;
         n    :: Int = 100,
         throw   :: Bool = false,
-        verbose :: Bool = false
+        verbose :: Bool = false,
+        sigdigits :: Int = 10
     ) -> tf :: Bool
 
 Check if a child GeoRegion defined by `cgeo` is within another GeoRegion `geo`.
@@ -80,6 +87,7 @@ Keyword Arguments
 - `n` : The number of segments to split each of the `GeoRegion`s into. Default is 100.
 - `throw`  : If `true`, then if `cgeo` is not within `geo`, an error is thrown and the program stops running.
 - `verbose` : If `true`, print logs to screen.
+- `sigdigits` : Specifies number of significant digits (i.e., precision) of the point coordinates used for checking. Defaults to 10.
 
 Returns
 =======
@@ -90,13 +98,20 @@ function Base.in(
     geo  :: GeoRegion;
     n    :: Int = 100,
     throw   :: Bool = false,
-    verbose :: Bool = false
+    verbose :: Bool = false,
+    sigdigits :: Int = 10
 )
 
     verbose ? (@info "$(modulelog()) - Performing a check to determine if the $(cgeo.name) GeoRegion ($(cgeo.ID)) is inside the $(geo.name) GeoRegion ($(geo.ID))") : nothing
 
     lon,lat = coordinates(cgeo,n=n)
-    isin = sum(.!in.(Point.(lon,lat),[geo]));
+
+    tlon,tlat = coordinates(geo,n=n)
+    tlon = round.(tlon,sigdigits=sigdigits)
+    tlat = round.(tlat,sigdigits=sigdigits)
+    tgeo = GeoRegion(tlon,tlat)
+
+    isin = sum(.!in.(Point.(lon,lat),[tgeo],sigdigits=sigdigits))
 
     if iszero(isin)
 

diff --git a/src/isin/ison.jl b/src/isin/ison.jl
@@ -2,7 +2,8 @@
     on(
         point :: Point2{<:Real},
         geo   :: GeoRegion;
-        throw :: Bool = false
+        throw :: Bool = false,
+        sigdigits :: Int = 10
     ) -> tf :: Bool
 
 Check if a geographical point `point` is on the boundary of a shape of a GeoRegion defined by `geo`.
@@ -15,6 +16,7 @@ Arguments
 Keyword Arguments
 =================
 - `throw` : If `true`, then if `point` is not within `geo`, an error is thrown and the program stops running.
+- `sigdigits` : Specifies number of significant digits (i.e., precision) of the point coordinates used for checking. Defaults to 10.
 
 Returns
 =======
@@ -23,21 +25,26 @@ Returns
 function on(
     point :: Point2{<:Real},
     geo   :: GeoRegion;
-    throw :: Bool = false
+    throw :: Bool = false,
+    sigdigits :: Int = 10
 )
 
     throw ? (@info "$(modulelog()) - Performing a check to determine if the coordinates $(point) are within the specified region boundaries.") : nothing
 
     plon = point[1]
-    plat = point[2]
+    plat = round(point[2],sigdigits=sigdigits)
 
     while plon > 360;  plon -= 360 end
     while plon < -180; plon += 360 end
 
+    plon1 = round(plon,sigdigits=sigdigits)
+    plon2 = round(plon-360,sigdigits=sigdigits)
+    plon3 = round(plon+360,sigdigits=sigdigits)
+
     isin = !iszero(sum([
-        touches(Point(plon    ,plat),geo.geometry.polygon),
-        touches(Point(plon+360,plat),geo.geometry.polygon),
-        touches(Point(plon-360,plat),geo.geometry.polygon)
+        touches(Point(plon1,plat),geo.geometry.polygon),
+        touches(Point(plon2,plat),geo.geometry.polygon),
+        touches(Point(plon3,plat),geo.geometry.polygon)
     ]))
 
     if !isin
@@ -62,7 +69,8 @@ end
         geo2 :: GeoRegion;
         n    :: Int = 2,
         throw   :: Bool = false,
-        verbose :: Bool = false
+        verbose :: Bool = false,
+        sigdigits :: Int = 10
     ) -> tf :: Bool
 
 Check if the GeoRegions `geo1` and `geo2` have the same shape. The order of `geo1` and `geo2` does not matter.
@@ -77,6 +85,7 @@ Keyword Arguments
 - `n` : The number of segments to split each of the `GeoRegion`s into. Default is 2.
 - `throw` : If `true`, then if `geo1` does not have the same shape as `geo2`, an error is thrown and the program stops running.
 - `verbose` : If `true`, print logs to screen.
+- `sigdigits` : Specifies number of significant digits (i.e., precision) of the point coordinates used for checking. Defaults to 10.
 
 Returns
 =======
@@ -87,14 +96,22 @@ function on(
     geo2 :: GeoRegion;
     n    :: Int = 2,
     throw   :: Bool = false,
-    verbose :: Bool = false
+    verbose :: Bool = false,
+    sigdigits :: Int = 10
 )
 
     verbose ? (@info "$(modulelog()) - Performing a check to determine if the $(geo1.name) GeoRegion  \"$(geo1.ID)\" shares the same shape as GeoRegion  \"$(geo2.ID)\".") : nothing
 
     lon1,lat1 = coordinates(geo1,n=n)
     lon2,lat2 = coordinates(geo2,n=n)
-    isin = sum(.!on.(Point.(lon1,lat1),[geo2])) + sum(.!on.(Point.(lon2,lat2),[geo1]))
+
+    tlon1 = round.(lon1,sigdigits=sigdigits); tlat1 = round.(lat1,sigdigits=sigdigits)
+    tlon2 = round.(lon2,sigdigits=sigdigits); tlat2 = round.(lat2,sigdigits=sigdigits)
+    tgeo1 = GeoRegion(tlon1,tlat1)
+    tgeo2 = GeoRegion(tlon2,tlat2)
+
+    isin = sum(.!on.(Point.(lon1,lat1),[tgeo2],sigdigits=sigdigits)) + 
+           sum(.!on.(Point.(lon2,lat2),[tgeo1],sigdigits=sigdigits))
 
     if iszero(isin)