Add polygon field to structs (#96)

* Swap cut polygons

* Fix removing cut poly

* Swap angle calculations

* Remove scaling polygons function

* Swap area function

* Remove calc_point_poly_dist

* Swap out points in poly

* Remove sort regions function

* Remove intersect_coords

* Rewrite moment of inertia calculations

* Finish separate_xy removal

* Switch basic calls to use get_polygons

* Add Extents

* Remove LG from grid output writer

* Debugging gridded data difference

* Remove LG from output

* Update poly_to_floes

* Finish removing LG clipping calls from all but floe utils

* Remove basic multipolygon functions

* Remove remaining LibGEOS intersection calls

* Remove bounding box utils

* Fix multipolygon clipping calls

* Abstract LG.Polygon calls out

* Abstract rest of LG functions  out

* Update ridge raft test

* Remove perturbation

* Small changes for testing

* Remove perturbation

* Readd perturb vec for random seed

* Remove unneeded print statements

* Remove LibGEOS from Subzero usage

* Add polygon field to struct

* Swap vertex finder code to use polygon

* Finish merge

* Update generate_subfloe_points to work with polygon

* Add poly field to topography and boundaries

* Clean up collision tests

* Use translate polygon helper function

* Remove unneeded make_polygon calls

* Remove rmholes function

* Update polygon rotation and movement

* Debug angles and optimize

* Add FT types to clipping

* Update grid cell polygon process

* Remove orient_coords function

* Remove remaining un-needed make_polygon calls

* Clean up Project.toml

* Small PR fixes

Project.toml

@@ -4,12 +4,11 @@ authors = ["Skylar Gering"]
 version = "0.1.0"
 
 [deps]
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
-GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 GeometryOps = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
@@ -21,19 +20,18 @@ NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 VoronoiCells = "e3e34ffb-84e9-5012-9490-92c94d0c60a4"
 
 [compat]
-BenchmarkTools = "1"
 CairoMakie = "0.12"
 DataStructures = "0.18"
 Extents = "0.1"
-GeometryBasics = "0.4"
-GeometryOps = "0.1.8"
+GeometryOps = "0.1.10"
 Interpolations = "0.14, 0.15"
 JLD2 = "0.4"
 Makie = "0.19, 0.20, 0.21"

examples/shear_flow.jl

@@ -44,7 +44,7 @@ floe_settings = FloeSettings(
 # Floe creation
 floe_arr = initialize_floe_field(
  FT,
- 50,
+ 500,
  [0.8],
  domain,
  hmean,
@@ -72,14 +72,13 @@ simulation = Simulation(
  model = model,
  consts = consts,
  Δt = Δt,
- nΔt = 2000,
+ nΔt = 5000,
  verbose = true,
  writers = writers,
  rng = Xoshiro(1),
  coupling_settings = coupling_settings,
 )
 run_time!(simulation)
-
 plot_sim(
  "output/shear_flow/floes.jld2",
  "output/shear_flow/initial_state.jld2",

examples/simple_strait.jl

@@ -87,7 +87,7 @@ simulation = Simulation(
  model = model,
  consts = consts,
  Δt = Δt,
- nΔt = 10000,
+ nΔt = 2500,
  verbose = true,
  floe_settings = floe_settings,
  coupling_settings = coupling_settings,

src/Subzero.jl

@@ -71,9 +71,11 @@ export
 import Base.@kwdef # this is being exported as of version 1.9
 import GeometryOps as GO
 import GeometryOps.GeoInterface as GI
-using CairoMakie, DataStructures, Dates, Extents, GeometryBasics, Interpolations, JLD2,
- LinearAlgebra, Logging, Makie, Measures, NCDatasets, NetCDF, Printf, Random,
- SplitApplyCombine, StaticArrays, Statistics, StructArrays, VoronoiCells
+import StaticArrays as SA
+using CairoMakie, CoordinateTransformations, DataStructures, Dates, Extents,
+ Interpolations, JLD2, LinearAlgebra, Logging, Makie, Measures, NCDatasets, NetCDF,
+ Printf, Random, Rotations, SplitApplyCombine, Statistics, StructArrays,
+ VoronoiCells
 
 """
 Coordinates are vector of vector of vector of points of the form:
@@ -99,6 +101,9 @@ const RingVec{T} = R where {
 
 const Polys{T} = GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, Nothing}}, Nothing, Nothing} where T <: AbstractFloat
 
+Base.convert(::Type{Polys{Float32}}, p::Polys{Float64}) = GO.tuples(p, Float32)
+Base.convert(::Type{Polys{Float64}}, p::Polys{Float32}) = GO.tuples(p, Float64)
+
 # Model
 include("simulation_components/floe.jl")
 include("floe_utils.jl")

src/floe_utils.jl

@@ -40,7 +40,7 @@ Input:
 Output:
  <PolyVec> representing the floe's coordinates xy plane
 """
-find_poly_coords(poly::Polys) = GI.coordinates(poly)
+find_poly_coords(poly) = GI.coordinates(poly)
 
 """
  intersect_polys(p1, p2)
@@ -52,19 +52,61 @@ Inputs:
 Output:
  Vector of Polygons
 """
-intersect_polys(p1, p2; kwargs...) = GO.intersection(p1, p2; target = GI.PolygonTrait(), fix_multipoly = nothing)
-diff_polys(p1, p2; kwargs...) = GO.difference(p1, p2; target = GI.PolygonTrait(), fix_multipoly = nothing) 
-union_polys(p1, p2; kwargs...) = GO.union(p1, p2; target = GI.PolygonTrait(), fix_multipoly = nothing)
-
+intersect_polys(p1, p2, ::Type{FT} = Float64; kwargs...) where FT = GO.intersection(p1, p2, FT; target = GI.PolygonTrait(), fix_multipoly = nothing)
+diff_polys(p1, p2, ::Type{FT} = Float64; kwargs...) where FT = GO.difference(p1, p2, FT; target = GI.PolygonTrait(), fix_multipoly = nothing) 
+union_polys(p1, p2, ::Type{FT} = Float64; kwargs...) where FT = GO.union(p1, p2, FT; target = GI.PolygonTrait(), fix_multipoly = nothing)
 simplify_poly(p, tol) = GO.simplify(p; tol = tol)
 
+function _translate_poly(::Type{FT}, p, Δx, Δy) where FT
+ t = CoordinateTransformations.Translation(Δx, Δy)
+ # TODO: can remove the tuples call after GO SVPoint PR
+ return GO.tuples(GO.transform(t, p), FT)
+end
+
+function _translate_floe!(::Type{FT}, floe, Δx, Δy) where FT
+ translate!(floe.coords, Δx, Δy)
+ floe.centroid[1] += Δx
+ floe.centroid[2] += Δy
+ floe.poly = _translate_poly(FT, floe.poly, Δx, Δy)
+ return
+end
+
+function _move_poly(::Type{FT}, poly, Δx, Δy, Δα, cx = zero(FT), cy = zero(FT)) where FT
+ rot = CoordinateTransformations.LinearMap(Rotations.Angle2d(Δα))
+ cent_rot = CoordinateTransformations.recenter(rot, (cx, cy))
+ trans = CoordinateTransformations.Translation(Δx, Δy)
+ # TODO: can remove the tuples call after GO SVPoint PR
+ return GO.tuples(GO.transform(trans ∘ cent_rot, poly), FT)
+end
+
+function _move_floe!(::Type{FT}, floe, Δx, Δy, Δα) where FT
+ cx, cy = floe.centroid
+ # Move coordinates and centroid
+ translate!(floe.coords, -cx, -cy)
+ rotate_radians!(floe.coords, Δα)
+ floe.centroid[1] += Δx
+ floe.centroid[2] += Δy
+ translate!(floe.coords, cx + Δx, cy + Δy)
+ # Move Polygon
+ floe.poly = _move_poly(FT, floe.poly, Δx, Δy, Δα, cx, cy)
+ return 
+end
+
 make_polygon(coords::PolyVec) = GI.Polygon(GO.tuples(coords))
 make_polygon(tuple_coords) = GI.Polygon(tuple_coords)
 make_polygon(ring::GI.LinearRing) = GI.Polygon([ring])
 make_multipolygon(coords::Vector{<:PolyVec}) = GI.MultiPolygon(GO.tuples(coords))
 make_multipolygon(tuple_coords) = GI.MultiPolygon(tuple_coords)
 make_multipolygon(polys::Vector{<:GI.Polygon}) = GI.MultiPolygon(polys)
 
+get_floe(floes::StructArray, i::Int) = LazyRow(floes, i)
+
+function _make_bounding_box_polygon(::Type{FT}, xmin, xmax, ymin, ymax) where FT
+ points = ((xmin, ymin), (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin))
+ ring = GI.LinearRing(SA.SVector{5, Tuple{FT, FT}}(points))
+ return GI.Polygon(SA.SVector(ring))
+end
+
 """
  deepcopy_floe(floe::LazyRow{Floe{FT}})
 
@@ -76,9 +118,11 @@ Outputs:
  they share values, but not referance.
 """
 function deepcopy_floe(floe::LazyRow{Floe{FT}}) where {FT}
+ poly = GO.tuples(floe.poly, FT)
  f = Floe{FT}(
+ poly = poly,
  centroid = copy(floe.centroid),
- coords = translate(floe.coords, 0, 0),
+ coords = find_poly_coords(poly),
  height = floe.height,
  area = floe.area,
  mass = floe.mass,
@@ -180,18 +224,6 @@ function rotate_radians!(coords::PolyVec, α)
  return
 end
 
-"""
- rotate_degrees!(coords::PolyVec, α)
-
-Rotate a polygon's coordinates by α degrees around the origin.
-Inputs:
- coords <PolyVec{AbstractFloat}> polygon coordinates
- α <Real> degrees to rotate the coordinates
-Outputs:
- Updates coordinates in place
-"""
-rotate_degrees!(coords::PolyVec, α) = rotate_radians!(coords, α * π/180)
-
 """
  hashole(coords::PolyVec{FT})
 
@@ -218,39 +250,16 @@ function hashole(poly::Polys)
  return GI.nhole(poly) > 0
 end 
 
-"""
- rmholes(coords::PolyVec{FT})
-
-Remove polygon coordinates's holes if they exist
-Inputs:
- coords <PolyVec{Float}> polygon coordinates
-Outputs:
- <PolyVec{Float}> polygonc coordinates after removing holes
-"""
-function rmholes(coords::PolyVec{FT}) where {FT<:AbstractFloat}
- return [coords[1]]
-end
 
 function rmholes!(coords::PolyVec{FT}) where {FT<:AbstractFloat}
  if length(coords) > 1
  deleteat!(coords, 2:length(coords))
  end
 end
 
-"""
- rmholes(poly::Polys)
-
-Remove polygon's holes if they exist
-Inputs:
- poly <Polygon> polygon
-Outputs:
- <Polygon> polygon without any holes
-"""
-function rmholes(poly::Polys)
- if hashole(poly)
- return make_polygon(GI.getexterior(poly))
- end
- return poly
+function rmholes!(poly::Polys)
+ deleteat!(poly.geom, 2:GI.nring(poly))
+ return
 end
 
 #=
@@ -293,7 +302,7 @@ function _calc_moment_inertia(
 end
 
 # Find the length of the maximum radius of a given polygon
-function _calc_max_radius(::Type{T}, poly, cent) where T
+function calc_max_radius(poly, cent, ::Type{T}) where T
  max_rad_sqrd = zero(T)
  Δx, Δy = GO._tuple_point(cent, T)
  for pt in GI.getpoint(GI.getexterior(poly))
@@ -307,110 +316,14 @@ function _calc_max_radius(::Type{T}, poly, cent) where T
  return sqrt(max_rad_sqrd)
 end
 
-"""
- orient_coords(coords)
-
-Take given coordinates and make it so that the first point has the smallest
-x-coordiante and so that the coordinates are ordered in a clockwise sequence.
-Duplicates vertices will be removed and the coordiantes will be closed (first
-and last point are the same).
-
-Input:
- coords <RingVec> vector of points [x, y]
-Output:
- coords <RingVec> oriented clockwise with smallest x-coordinate first
-"""
-function orient_coords(coords::RingVec)
- # extreem_idx is point with smallest x-value - if tie, choose lowest y-value
- extreem_idx = 1
- for i in eachindex(coords)
- ipoint = coords[i]
- epoint = coords[extreem_idx]
- if ipoint[1] < epoint[1]
- extreem_idx = i
- elseif ipoint[1] == epoint[1] && ipoint[2] < epoint[2]
- extreem_idx = i
- end
- end
- # extreem point must be first point in list
- new_coords = similar(coords)
- circshift!(new_coords, coords, -extreem_idx + 1)
- valid_ringvec!(new_coords)
-
- # if coords are counterclockwise, switch to clockwise
- orient_matrix = hcat(
- ones(3),
- vcat(new_coords[1]', new_coords[2]', new_coords[end-1]') # extreem/adjacent points
- )
- if det(orient_matrix) > 0
- reverse!(new_coords)
- end
- return new_coords
-end
-
-"""
- intersect_lines(l1, l2)
-
-Finds the intersection points of two curves l1 and l2. The curves l1, l2 can be
-either closed or open. In this version, l1 and l2 must be distinct. If no
-intersections are found, the returned P is empty.
-Inputs:
- l1 <PolyVec{Float}> line/polygon coordinates
- l2 <PolyVec{Float}> line/polygon coordinates
-Outputs:
- <Set{Tuple{Float, Float}}> Set of points that are at the intersection of the
- two line segments.
-"""
-function intersect_lines(l1::PolyVec{FT}, l2) where {FT}
- points = Vector{Tuple{FT, FT}}()
- for i in 1:(length(l1[1]) - 1)
- # First line represented as a1x + b1y = c1
- x11 = l1[1][i][1]
- x12 = l1[1][i+1][1]
- y11 = l1[1][i][2]
- y12 = l1[1][i+1][2]
- a1 = y12 - y11
- b1 = x11 - x12
- c1 = a1 * x11 + b1 * y11
- for j in 1:(length(l2[1]) - 1)
- # Second line represented as a2x + b2y = c2
- x21 = l2[1][j][1]
- x22 = l2[1][j+1][1]
- y21 = l2[1][j][2]
- y22 = l2[1][j+1][2]
- a2 = y22 - y21
- b2 = x21 - x22
- c2 = a2 * x21 + b2 * y21
-
- determinant = a1 * b2 - a2 * b1
- # Find place there two lines cross
- # Note that lines extend beyond given line segments
- if determinant != 0
- x = (b2*c1 - b1*c2)/determinant
- y = (a1*c2 - a2*c1)/determinant
- p = (x, y)
- # Make sure intersection is on given line segments
- if min(x11, x12) <= x <= max(x11, x12) &&
- min(x21, x22) <= x <= max(x21, x22) &&
- min(y11, y12) <= y <= max(y11, y12) &&
- min(y21, y22) <= y <= max(y21, y22) &&
- !(p in points)
- push!(points, p)
- end
- end
- end
- end
- return points
-end
-
 """
  which_vertices_match_points(ipoints, coords, atol)
 
 Find which vertices in coords match given points
 Inputs:
  points <Vector{Tuple{Float, Float} or Vector{Vector{Float}}}> points to
  match to vertices within polygon
- coords <PolVec> polygon coordinates
+ region <Polygon> polygon 
  atol <Float> distance vertex can be away from target point before being
  classified as different points
 Output:
@@ -419,24 +332,17 @@ Note:
  If last coordinate is a repeat of first coordinate, last coordinate index is
  NOT recorded.
 """
-function which_vertices_match_points(
- points,
- coords::PolyVec{FT},
- atol = 1,
-) where {FT}
+function which_vertices_match_points(points, region::Polys{FT}, atol = 1) where FT
  idxs = Vector{Int}()
  npoints = length(points)
  if points[1] == points[end]
  npoints -= 1
  end
- @views for i in 1:npoints # find which vertex matches point
+ for i in 1:npoints # find which vertex matches point
  min_dist = FT(Inf)
  min_vert = 1
- for j in eachindex(coords[1])
- dist = sqrt(
- (coords[1][j][1] - points[i][1])^2 +
- (coords[1][j][2] - points[i][2])^2,
- )
+ for (j, pt) in enumerate(GI.getpoint(GI.getexterior(region)))
+ dist = sqrt(GO.distance(pt, points[i], FT))
  if dist < min_dist
  min_dist = dist
  min_vert = j