Merge pull request #193 from gaelforget/polygons_202512

Polygons operations

Author: gaelforget

Project.toml

@@ -1,6 +1,6 @@
 name = "MeshArrays"
 uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
-version = "0.4.0"
+version = "0.4.1"
 authors = ["gaelforget <[email protected]>"]
 
 [deps]
@@ -21,6 +21,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [weakdeps]
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
+GeometryOps = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
@@ -29,6 +30,7 @@ Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 [extensions]
 MeshArraysDataDepsExt = ["DataDeps"]
 MeshArraysGeoJSONExt = ["GeoJSON"]
+MeshArraysGeometryOpsExt = ["GeometryOps"]
 MeshArraysJLD2Ext = ["JLD2"]
 MeshArraysMakieExt = ["Makie"]
 MeshArraysProjExt = ["Proj"]
@@ -39,6 +41,7 @@ CatViews = "1.0"
 DataDeps = "0.7"
 GeoInterface = "1.4, 1.5, 1.6"
 GeoJSON = "0.6, 0.7, 0.8"
+GeometryOps = "0.1"
 Glob = "1"
 JLD2 = "0.4, 0.5, 0.6"
 Makie = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24"

examples/basics.jl

@@ -94,7 +94,7 @@ end
 
 # ╔═╡ 4a3f3b75-1c75-4ce5-8e75-f4dd263199cd
 let
-	γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
+	γ=GridSpec(ID=:LLC90())
 
 	a=fill(2.0,MeshArray(γ))
 	aa=write(a)
@@ -158,15 +158,15 @@ md"""## 2. Cube Sphere Grid
 
 Now, let's instead use a grid that has 6 subdomains, 32x32 points each, covering the six faces of a cube. This _cube sphere_ grid has connections between subdomains that are slightly more complicated than in the first example. 
 
-Here, grid variables are read from the `MeshArrays.GRID_CS32` folder by the `GridLoad` function. The `GridSpec` function provides the corresponding subdomain sizes. Again we call `comp()` which combines steps 2 and 3, and `viz()` to visualize. 
+Here, grid variables are read from the `GRID_CS32` folder by the `GridLoad` function. The `GridSpec` function provides the corresponding subdomain sizes. Again we call `comp()` which combines steps 2 and 3, and `viz()` to visualize. 
 """
 
 # ╔═╡ d37ae153-8646-44ae-abd9-7e0e7d479275
 md"""## 3. LLC90 grid
 
 The [Lat-Lon-Cap grid](http://www.geosci-model-dev.net/8/3071/2015/) (or LLC) is a global ocean model grid which is widely used in the [MITgcm user community](https://mitgcm.readthedocs.io/en/latest/). It has 5 uneven subdomains, variable grid spacing, and continents [(Forget et al 2015)](http://www.geosci-model-dev.net/8/3071/2015/). LLC90's resolution is one degree albeit with modications in the Arctic and along the Equator.
 
-Here, grid variables are read from the `MeshArrays.GRID_LLC90` folder by the `GridLoad` function. The `GridSpec` function provides the corresponding subdomain sizes. Again we call `comp()` which combines steps 2 and 3, and `viz()` to visualize. 
+Here, grid variables are read from the `GRID_LLC90` folder by the `GridLoad` function. The `GridSpec` function provides the corresponding subdomain sizes. Again we call `comp()` which combines steps 2 and 3, and `viz()` to visualize. 
 """
 
 # ╔═╡ 0865dd0f-c43e-43a7-aa65-74aba4f4460d
@@ -248,7 +248,7 @@ end
 
 # ╔═╡ dacc2261-1d75-4d68-b82c-174e4fae3631
 begin
-	γ_b=GridSpec("CubeSphere",MeshArrays.GRID_CS32)
+	γ_b=GridSpec("CubeSphere",MeshArrays.Dataset("GRID_CS32"))
 	Γ_b=GridLoad(γ_b;option="full")
 
 	Rini_b,Rend_b=smoothing_example(Γ_b)
@@ -265,7 +265,7 @@ end
 
 # ╔═╡ 23198b8d-6bf0-4af9-91d0-7886c9574f81
 begin
-	γ_c=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
+	γ_c=GridSpec(ID=:LLC90)
 	Γ_c=GridLoad(γ_c;option="full")
 	Rini_c,Rend_c=smoothing_example(Γ_c)
 	PlutoUI.with_terminal() do
@@ -777,18 +777,36 @@ git-tree-sha1 = "7a214fdac5ed5f59a22c2d9a885a16da1c74bbc7"
 uuid = "559328eb-81f9-559d-9380-de523a88c83c"
 version = "1.0.17+0"
 
+[[deps.GeoFormatTypes]]
+git-tree-sha1 = "7528a7956248c723d01a0a9b0447bf254bf4da52"
+uuid = "68eda718-8dee-11e9-39e7-89f7f65f511f"
+version = "0.4.5"
+
+[[deps.GeoInterface]]
+deps = ["DataAPI", "Extents", "GeoFormatTypes"]
+git-tree-sha1 = "b7c5cdf45298877bb683bdda3f871ff7070985c4"
+uuid = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
+version = "1.6.0"
+
+    [deps.GeoInterface.extensions]
+    GeoInterfaceMakieExt = ["Makie", "GeometryBasics"]
+    GeoInterfaceRecipesBaseExt = "RecipesBase"
+
+    [deps.GeoInterface.weakdeps]
+    GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+    Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+    RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+
 [[deps.GeometryBasics]]
 deps = ["EarCut_jll", "Extents", "IterTools", "LinearAlgebra", "PrecompileTools", "Random", "StaticArrays"]
 git-tree-sha1 = "1f5a80f4ed9f5a4aada88fc2db456e637676414b"
 uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 version = "0.5.10"
+weakdeps = ["GeoInterface"]
 
     [deps.GeometryBasics.extensions]
     GeometryBasicsGeoInterfaceExt = "GeoInterface"
 
-    [deps.GeometryBasics.weakdeps]
-    GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
-
 [[deps.GettextRuntime_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll"]
 git-tree-sha1 = "45288942190db7c5f760f59c04495064eedf9340"
@@ -1221,10 +1239,10 @@ uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
 version = "0.6.7"
 
 [[deps.MeshArrays]]
-deps = ["CatViews", "Dates", "Distributed", "Glob", "LazyArtifacts", "NearestNeighbors", "Pkg", "Printf", "SharedArrays", "SparseArrays", "Statistics", "Unitful"]
-git-tree-sha1 = "3ea2dc9aaaa4f4aa8799d568a4fef5abcfd2b7bf"
+deps = ["CatViews", "Dates", "Distributed", "GeoInterface", "Glob", "LazyArtifacts", "NearestNeighbors", "Pkg", "Printf", "SharedArrays", "SparseArrays", "Statistics", "Unitful"]
+git-tree-sha1 = "bc4820b8c8648fc0d60503562a6acae0aeef312d"
 uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
-version = "0.3.24"
+version = "0.4.0"
 
     [deps.MeshArrays.extensions]
     MeshArraysDataDepsExt = ["DataDeps"]

examples/geography.jl

@@ -125,7 +125,7 @@ md"""## Appendices
 
 # ╔═╡ 39924391-38ce-46a1-877f-80a7975340a0
 begin
-	pth=MeshArrays.GRID_LLC90
+    pth=MeshArrays.Dataset("GRID_LLC90")
 	γ=GridSpec("LatLonCap",pth)
 	Γ=GridLoad(γ;option="light")	
 	#LC=LatitudeCircles(-89.0:89.0,Γ)
@@ -248,9 +248,9 @@ end
 # ╔═╡ 229d395f-f3b4-40df-a482-264d540be875
 begin
 	if source==:shp_example
-		fil=demo.download_polygons("ne_110m_admin_0_countries.shp")
+		fil=MeshArrays.Dataset("countries_shp1",do_read=false)
 	elseif source==:json_example
-		fil=demo.download_polygons("countries.geojson")
+		fil=MeshArrays.Dataset("countries_geojson1",do_read=false)
 	end
 	pol=MeshArrays.read_polygons(fil)
 	"Done With Reading Country Polygons"
@@ -260,7 +260,7 @@ end
 begin
 	fig_polygons=Figure()
 	ax1=Axis(fig_polygons[1,1])
-	[lines!(ax1,l1,color = :black, linewidth = 0.5) for l1 in pol]
+	[lines!(ax1,l1.geometry,color = :black, linewidth = 0.5) for l1 in pol]
 	fig_polygons
 end
 
@@ -506,9 +506,9 @@ version = "1.0.0"
 
 [[deps.Climatology]]
 deps = ["CSV", "DataDeps", "DataFrames", "DataStructures", "Dataverse", "Dates", "Distributed", "Glob", "JLD2", "MeshArrays", "Pkg", "Printf", "RollingFunctions", "STAC", "Scratch", "SharedArrays", "Statistics", "TOML", "URIs"]
-git-tree-sha1 = "12e7b88c1154042f725d914886604e92d59ba29a"
+git-tree-sha1 = "6919ec5589a0424cb1c81b40966f6ae9f0e89ac5"
 uuid = "9e9a4d37-2d2e-41e3-8b85-f7978328d9c7"
-version = "0.5.16"
+version = "0.5.17"
 
     [deps.Climatology.extensions]
     ClimatologyMITgcmExt = ["MITgcm"]
@@ -928,9 +928,9 @@ version = "5.2.3+0"
 
 [[deps.Glib_jll]]
 deps = ["Artifacts", "GettextRuntime_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "50c11ffab2a3d50192a228c313f05b5b5dc5acb2"
+git-tree-sha1 = "6b4d2dc81736fe3980ff0e8879a9fc7c33c44ddf"
 uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
-version = "2.86.0+0"
+version = "2.86.2+0"
 
 [[deps.Glob]]
 git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496"
@@ -1425,10 +1425,10 @@ uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
 version = "2.28.1010+0"
 
 [[deps.MeshArrays]]
-deps = ["CatViews", "Dates", "Distributed", "Glob", "LazyArtifacts", "NearestNeighbors", "Pkg", "Printf", "SharedArrays", "SparseArrays", "Statistics", "Unitful"]
-git-tree-sha1 = "3ea2dc9aaaa4f4aa8799d568a4fef5abcfd2b7bf"
+deps = ["CatViews", "Dates", "Distributed", "GeoInterface", "Glob", "LazyArtifacts", "NearestNeighbors", "Pkg", "Printf", "SharedArrays", "SparseArrays", "Statistics", "Unitful"]
+git-tree-sha1 = "bc4820b8c8648fc0d60503562a6acae0aeef312d"
 uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
-version = "0.3.24"
+version = "0.4.0"
 weakdeps = ["DataDeps", "GeoJSON", "JLD2", "Makie", "Proj", "Shapefile"]
 
     [deps.MeshArrays.extensions]

examples/vectors.jl

@@ -84,7 +84,7 @@ Here we load `MeshArrays.jl` and packages that provide a climatological estimate
 
 # ╔═╡ 3c2de9d4-b091-464b-9210-aa84f3d4c5f1
 begin
-	pth=MeshArrays.GRID_LLC90
+	pth=MeshArrays.Dataset("GRID_LLC90")
 	γ=GridSpec("LatLonCap",pth)
 	Γ=GridLoad(γ;option="full")
 	LC=LatitudeCircles(-89.0:89.0,Γ)
@@ -521,9 +521,9 @@ version = "0.3.11"
 
 [[deps.Climatology]]
 deps = ["CSV", "DataDeps", "DataFrames", "DataStructures", "Dataverse", "Dates", "Distributed", "Glob", "JLD2", "MeshArrays", "Pkg", "Printf", "RollingFunctions", "STAC", "Scratch", "SharedArrays", "Statistics", "TOML", "URIs"]
-git-tree-sha1 = "12e7b88c1154042f725d914886604e92d59ba29a"
+git-tree-sha1 = "6919ec5589a0424cb1c81b40966f6ae9f0e89ac5"
 uuid = "9e9a4d37-2d2e-41e3-8b85-f7978328d9c7"
-version = "0.5.16"
+version = "0.5.17"
 
     [deps.Climatology.extensions]
     ClimatologyMITgcmExt = ["MITgcm"]
@@ -960,9 +960,9 @@ version = "2.52.0+0"
 
 [[deps.Glib_jll]]
 deps = ["Artifacts", "GettextRuntime_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "50c11ffab2a3d50192a228c313f05b5b5dc5acb2"
+git-tree-sha1 = "6b4d2dc81736fe3980ff0e8879a9fc7c33c44ddf"
 uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
-version = "2.86.0+0"
+version = "2.86.2+0"
 
 [[deps.Glob]]
 git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496"
@@ -1425,9 +1425,9 @@ version = "1.1.0"
 
 [[deps.MITgcm]]
 deps = ["ClimateModels", "DataDeps", "Dataverse", "Dates", "Distributed", "FortranFiles", "Glob", "MeshArrays", "Printf", "Scratch", "SharedArrays", "SparseArrays", "Statistics", "StyledStrings", "UUIDs"]
-git-tree-sha1 = "38eaef144504c69593362406b4425d3d2d1a0c9a"
+git-tree-sha1 = "abce854796b41317d2a33eea8714b27cb658cc14"
 uuid = "dce5fa8e-68ce-4431-a242-9469c69627a0"
-version = "0.5.11"
+version = "0.5.12"
 weakdeps = ["NetCDF"]
 
     [deps.MITgcm.extensions]
@@ -1503,10 +1503,10 @@ uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
 version = "2.28.1010+0"
 
 [[deps.MeshArrays]]
-deps = ["CatViews", "Dates", "Distributed", "Glob", "LazyArtifacts", "NearestNeighbors", "Pkg", "Printf", "SharedArrays", "SparseArrays", "Statistics", "Unitful"]
-git-tree-sha1 = "3ea2dc9aaaa4f4aa8799d568a4fef5abcfd2b7bf"
+deps = ["CatViews", "Dates", "Distributed", "GeoInterface", "Glob", "LazyArtifacts", "NearestNeighbors", "Pkg", "Printf", "SharedArrays", "SparseArrays", "Statistics", "Unitful"]
+git-tree-sha1 = "bc4820b8c8648fc0d60503562a6acae0aeef312d"
 uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
-version = "0.3.24"
+version = "0.4.0"
 
     [deps.MeshArrays.extensions]
     MeshArraysDataDepsExt = ["DataDeps"]

ext/MeshArraysGeometryOpsExt.jl

@@ -0,0 +1,36 @@
+
+############################################################
+#                                                          #
+#   Geometric Operations on Polygons                       #
+#                                                          #
+############################################################
+
+module MeshArraysGeometryOpsExt
+
+    import GeometryOps as GO
+    import MeshArrays: within_pol, GI
+
+    """
+        within_pol(pol; ID=1)
+
+    Generate a `name,rule` pair to test if location 
+        `lon,lat` is within polygon `pol[ID].geometry`.
+
+    ```
+    using MeshArrays, DataDeps, GeoJSON, GeometryOps
+    pol=MeshArrays.Dataset("oceans_geojson1")
+    name,rule=MeshArrays.within_pol(pol; ID=11)
+    name,rule(-30,40)
+    ```
+    """
+    function within_pol(pol; ID=0)
+        name,geom=if ID==0
+            "unknown",pol
+        else
+            pol[ID].name,pol[ID].geometry
+        end
+        rule = (x,y) -> GO.within(GI.Point(x,y), geom)
+        (name,rule)
+    end
+
+end #module MeshArraysGeometryOpsExt

src/MeshArrays.jl

@@ -14,6 +14,7 @@ function read_jld2 end
 function write_jld2 end
 function ProjAxis end
 function grid_lines! end
+function within_pol end
 
 include("Types.jl")
 include("Grids.jl")

test/Project.toml

@@ -3,6 +3,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
+GeometryOps = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"

test/runtests.jl

@@ -1,13 +1,13 @@
 using Test, Documenter, Suppressor, MeshArrays, CairoMakie
-import DataDeps, JLD2, Shapefile, GeoJSON, Proj
+import DataDeps, JLD2, Shapefile, GeoJSON, Proj, GeometryOps
 
 MeshArrays.Dataset("GRID_LL360")
 MeshArrays.Dataset("GRID_LLC90")
 MeshArrays.Dataset("GRID_LLC270")
 MeshArrays.Dataset("GRID_CS32")
 
 pol_shp=MeshArrays.Dataset("countries_shp1")
-pol_json=MeshArrays.Dataset("countries_geojson1")
+pol_json=MeshArrays.Dataset("oceans_geojson1")
 
 p=dirname(pathof(MeshArrays))
 include(joinpath(p,"../examples/Demos.jl"))
@@ -349,6 +349,14 @@ end
     plot_examples(:baseproj,proj,lon0,pol=pol_shp)
 end
 
+@testset "polygon operattions" begin
+    name,rule=MeshArrays.within_pol(pol_json; ID=11)
+    rule_vec = (x,y) -> rule.(x,y)
+
+    np=10000; lo=-180 .+360*rand(np); la=-90 .+180*rand(np);
+    sum(rule_vec(lo,la))
+end
+
 @testset "doctests" begin
     doctest(MeshArrays; manual = false)
 end