Adding isin a TiltRegion

GeoRegionsEcosystem · Jul 26, 2024 · c5f4e70 · c5f4e70
1 parent 2905def
commit c5f4e70
Showing 1 changed file with 100 additions and 72 deletions.
diff --git a/src/isin/georegion.jl b/src/isin/georegion.jl
@@ -1,136 +1,157 @@
 """
     in(
         Child  :: GeoRegion,
-        polyG  :: PolyRegion;
-        domask :: Bool = false,
+        rectG  :: RectRegion;
         throw  :: Bool = false
     ) -> Bool
 
-Check if a child GeoRegion defined by `Child` is within a PolyRegion `polyG`.
+Check if a child GeoRegion defined by `Child` is within a RectRegion `rectG`.
 
 Arguments
 =========
 
 - `Child` : A GeoRegion that we postulate to be a "child", or a subset of the GeoRegion defined by `polyG`
-- `polyG` : A GeoRegion that we postulate to be a "parent", or containing the GeoRegion defined by `Child`
+- `rectG` : A GeoRegion that we postulate to be a "parent", or containing the GeoRegion defined by `Child`
 
 Keyword Arguments
 =================
 
 - `throw`  : If `true`, then if `Child` is not within `polyG`, an error is thrown and the program stops running
-- `domask` : If `throw` is `false` and `domask` is `true`, return a mask (with bounds defined by the `Child` GeoRegion) showing the region where `Child` and `polyG` do not overlap
 """
 Base.in(
-    Child  :: GeoRegion,
-    polyG  :: PolyRegion;
-    domask :: Bool = false,
-    throw  :: Bool = false
-) = isinGeoRegion(Child,polyG,domask=domask,throw=throw)
+    Child :: GeoRegion,
+    rectG :: RectRegion;
+    throw :: Bool = false
+) = isinGeoRegion(Child,rectG,throw=throw)
 
 """
     in(
-        Child  :: GeoRegion,
-        rectG  :: RectRegion;
-        throw  :: Bool = false
+        cgeo  :: GeoRegion,
+        geo   :: Union{TiltRegion,PolyRegion};
+        n     :: Int = 100,
+        throw :: Bool = false
     ) -> Bool
 
-Check if a child GeoRegion defined by `Child` is within a RectRegion `rectG`.
+Check if a child GeoRegion defined by `cgeo` is within a TiltRegion or PolyRegion `geo`.
 
 Arguments
 =========
 
-- `Child` : A GeoRegion that we postulate to be a "child", or a subset of the GeoRegion defined by `polyG`
-- `polyG` : A GeoRegion that we postulate to be a "parent", or containing the GeoRegion defined by `Child`
+- `cgeo` : A GeoRegion that we postulate to be a "child", or a subset of the GeoRegion defined by `polyG`
+- `geo` : A TiltRegion or PolyRegion that we postulate to be a "parent", or containing the GeoRegion defined by `Child`
 
 Keyword Arguments
 =================
 
-- `throw`  : If `true`, then if `Child` is not within `polyG`, an error is thrown and the program stops running
+- `n`  : Number of points per polygon side to test
+- `throw` : If `true`, then if `cgeo` is not within `geo`, an error is thrown and the program stops running
 """
 Base.in(
-    Child :: GeoRegion,
-    rectG :: RectRegion;
+    cgeo  :: GeoRegion,
+    geo   :: Union{TiltRegion,PolyRegion};
+    n     :: Int = 100,
     throw :: Bool = false
-) = isinGeoRegion(Child,rectG,throw=throw)
+) = isinGeoRegion(cgeo,geo,n=n,throw=throw)
 
 """
     isinGeoRegion(
         Child  :: GeoRegion,
-        polyG  :: PolyRegion;
-        domask :: Bool = false,
+        rectG  :: RectRegion;
         throw  :: Bool = true
     ) -> Bool
 
-Check if a child GeoRegion defined by `Child` is within a PolyRegion `polyG`.
+Check if a child GeoRegion defined by `Child` is within a RectRegion `rectG`.
 
 Arguments
 =========
 
 - `Child` : A GeoRegion that we postulate to be a "child", or a subset of the GeoRegion defined by `polyG`
-- `polyG` : A GeoRegion that we postulate to be a "parent", or containing the GeoRegion defined by `Child`
+- `rectG` : A GeoRegion that we postulate to be a "parent", or containing the GeoRegion defined by `Child`
 
 Keyword Arguments
 =================
 
-- `throw`  : If `true`, then if `Child` is not within `polyG`, an error is thrown and the program stops running
-- `domask` : If `throw` is `false` and `domask` is `true`, return a mask (with bounds defined by the `Child` GeoRegion) showing the region where `Child` and `polyG` do not overlap
+- `throw`  : If `true`, then if `Child` is not within `rectG`, an error is thrown and the program stops running
 """
 function isinGeoRegion(
-    Child  :: GeoRegion,
-    polyG  :: PolyRegion;
-    domask :: Bool = false,
-    throw  :: Bool = true
+    Child :: GeoRegion,
+    rectG :: RectRegion;
+    throw :: Bool = true,
 )
 
-    @info "$(modulelog()) - Performing a check to determine if the $(Child.name) GeoRegion ($(Child.ID)) is inside the $(polyG.name) GeoRegion ($(polyG.ID))"
-
-    N = Child.N
-    S = Child.S
-    E = Child.E
-    W = Child.W
-
-    lon = range(W,E,length=10001); nlon = length(lon)
-    lat = range(S,N,length=10001); nlat = length(lat)
-
-    mask = zeros(Bool,nlon,nlat)
+    @info "$(modulelog()) - Performing a check to determine if the $(Child.name) GeoRegion ($(Child.ID)) is inside the $(rectG.name) GeoRegion ($(rectG.ID))"
 
-    for ilat in 1 : nlat, ilon = 1 : nlon
-        if in(Point2(lon[ilon],lat[ilat]),Child)
-            if in(Point2(lon[ilon],lat[ilat]),polyG)
-                mask[ilon,ilat] = 1
-            end
-        end
-    end
+    isin = isgridinregion(
+        [Child.N,Child.S,Child.E,Child.W],
+        [rectG.N,rectG.S,rectG.E,rectG.W],
+        throw=false
+    )
 
-    if sum(mask) > 0
+    if !isin
 
         if throw
-            error("$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is not a subset of the GeoRegion $(polyG.ID) ($(polyG.name))")
+            error("$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is not a subset of the GeoRegion $(rectG.ID) ($(rectG.name))")
         else
-            if domask
-                @warn "$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is not a subset of the GeoRegion $(polyG.ID) ($(polyG.name)), returning a mask to show which regions do not intersect"
-                return mask
-            else
-                @warn "$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is not a subset of the GeoRegion $(polyG.ID) ($(polyG.name))"
-                return false
-            end
+            @warn "$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is not a subset of the GeoRegion $(rectG.ID) ($(rectG.name))"
+            return false
         end
 
     else
-        @info "$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is indeed a subset of the GeoRegion $(polyG.ID) ($(polyG.name))"
+        @info "$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is indeed a subset of the GeoRegion $(rectG.ID) ($(rectG.name))"
         return true
     end
 
 end
 
+
+
+"""
+    isinGeoRegion(
+        Child :: GeoRegion,
+        tiltG :: TiltRegion;
+        n     :: Int = 100,
+        throw :: Bool = false
+    ) -> Bool
+
+Check if a child GeoRegion defined by `Child` is within a TiltRegion `tiltG`.
+
+Arguments
+=========
+
+- `Child` : A GeoRegion that we postulate to be a "child", or a subset of the GeoRegion defined by `polyG`
+- `tiltG` : A GeoRegion that we postulate to be a "parent", or containing the GeoRegion defined by `Child`
+
+Keyword Arguments
+=================
+
+- `n`  : Number of points per polygon side to test
+- `throw`  : If `true`, then if `Child` is not within `tiltG`, an error is thrown and the program stops running
+"""
+function isinGeoRegion(
+    Child :: GeoRegion,
+    tiltG :: TiltRegion;
+    n     :: Int = 100,
+    throw :: Bool = true,
+)
+
+    @info "$(modulelog()) - Performing a check to determine if the $(Child.name) GeoRegion ($(Child.ID)) is inside the $(rectG.name) GeoRegion ($(rectG.ID))"
+
+    lon,lat = getTiltShape(tiltG)
+    tgeo = PolyRegion("","GLB","",lon,lat,save=false,verbose=false)
+
+    isinGeoRegion(Child,tgeo,n=n,throw=throw)
+
+end
+
 """
     isinGeoRegion(
         Child  :: GeoRegion,
-        rectG  :: RectRegion;
+        polyG  :: PolyRegion;
+        domask :: Bool = false,
         throw  :: Bool = true
     ) -> Bool
 
-Check if a child GeoRegion defined by `Child` is within a RectRegion `rectG`.
+Check if a child GeoRegion defined by `Child` is within a PolyRegion `polyG`.
 
 Arguments
 =========
@@ -141,33 +162,40 @@ Arguments
 Keyword Arguments
 =================
 
+- `n`  : Number of points per polygon side to test
 - `throw`  : If `true`, then if `Child` is not within `polyG`, an error is thrown and the program stops running
 """
 function isinGeoRegion(
     Child :: GeoRegion,
-    rectG :: RectRegion;
-    throw :: Bool = true
+    polyG :: PolyRegion;
+    n     :: Int = 100,
+    throw :: Bool = true,
 )
 
-    @info "$(modulelog()) - Performing a check to determine if the $(Child.name) GeoRegion ($(Child.ID)) is inside the $(rectG.name) GeoRegion ($(rectG.ID))"
+    @info "$(modulelog()) - Performing a check to determine if the $(Child.name) GeoRegion ($(Child.ID)) is inside the $(polyG.name) GeoRegion ($(polyG.ID))"
 
-    isin = isgridinregion(
-        [Child.N,Child.S,Child.E,Child.W],
-        [rectG.N,rectG.S,rectG.E,rectG.W],
-        throw=false
-    )
+    lon,lat = coordGeoRegion(Child,n=n)
+    npts = length(lon)
+
+    isin = 0
 
-    if !isin
+    for ipnt in npts
+        if in(Point2(lon[ipnt],lat[ipnt]),polyG)
+            isin += 1
+        end
+    end
+
+    if isin > 0
 
         if throw
-            error("$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is not a subset of the GeoRegion $(rectG.ID) ($(rectG.name))")
+            error("$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is not a subset of the GeoRegion $(polyG.ID) ($(polyG.name))")
         else
-            @warn "$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is not a subset of the GeoRegion $(rectG.ID) ($(rectG.name))"
+            @warn "$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is not a subset of the GeoRegion $(polyG.ID) ($(polyG.name))"
             return false
         end
 
     else
-        @info "$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is indeed a subset of the GeoRegion $(rectG.ID) ($(rectG.name))"
+        @info "$(modulelog()) - The GeoRegion $(Child.ID) ($(Child.name)) is indeed a subset of the GeoRegion $(polyG.ID) ($(polyG.name))"
         return true
     end