Feat(tracker): add rotation calculations

IFTPipeline.jl/Project.toml

@@ -7,6 +7,8 @@ version = "0.1.0"
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CSVFiles = "5d742f6a-9f54-50ce-8119-2520741973ca"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 Folds = "41a02a25-b8f0-4f67-bc48-60067656b558"
@@ -20,13 +22,17 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"
 
 [compat]
 ArgParse = "1.1"
+DataFrames = "1.7.0"
+Dates = "1.10.9"
 Folds = "0.2"
 HDF5 = "0.16, 0.17"
 IceFloeTracker = "0.6.4"
 LoggingExtras = "1.0"
 Pkg = "1.9"
 PyCall = "1.96"
 TOML = "1.0"
+TimeZones = "1.21.3"

IFTPipeline.jl/src/IFTPipeline.jl

@@ -49,6 +49,7 @@ include("preprocess.jl")
 include("feature-extraction.jl")
 include("tracker.jl")
 include("h5.jl")
+include("rotations.jl")
 
 export cache_vector,
     sharpen,
@@ -77,7 +78,8 @@ export cache_vector,
     mkclitrack!,
     mkfilenames,
     makeh5files,
-    getlatlon
+    getlatlon,
+    get_rotation_single
 
 export IceFloeTracker
 end

IFTPipeline.jl/src/cli.jl

@@ -4,6 +4,7 @@ using ArgParse
 using IceFloeTracker
 using IFTPipeline
 using Serialization
+using TimeZones
 
 function mkclipreprocess!(settings)
     @add_arg_table! settings["preprocess"] begin
@@ -318,7 +319,7 @@ function mkclimakeh5_single!(settings)
         "--passtime"
         help = "Satellite pass time"
         required = true
-        arg_type = DateTime
+        arg_type = ZonedDateTime
 
         "--truecolor"
         help = "Path to truecolor image"
@@ -498,7 +499,7 @@ function mkclitrack_single!(settings)
         help = "Path to object with satellite pass times"
         required = true
         nargs = '+'
-        arg_type = DateTime
+        arg_type = ZonedDateTime
 
         "--latlon"
         help = "Path to geotiff image with latitude/longitude data"
@@ -645,6 +646,19 @@ function mkclilandmask_single!(settings)
     return nothing
 end
 
+function mkcli_get_rotation_single!(settings)
+    @add_arg_table! settings["get_rotation_single"] begin
+        "--input", "-i"
+        help = "Tracked floes CSV file with a 'satellite' column"
+        required = true
+
+        "--output", "-o"
+        help = "Tracked floes CSV file with rotations calculated"
+        required = true
+    end
+    return nothing
+end
+
 function mkcli!(settings, common_args)
     d = Dict(
         "landmask" => mkclilandmask!,
@@ -658,6 +672,7 @@ function mkcli!(settings, common_args)
         "makeh5files_single" => mkclimakeh5_single!,
         "track" => mkclitrack!,
         "track_single" => mkclitrack_single!,
+        "get_rotation_single" => mkcli_get_rotation_single!,
     )
 
     for t in keys(d)
@@ -714,6 +729,10 @@ function main()
         "track_single"
         help = "Pair ice floes in day k with ice floes in day k+1"
         action = :command
+
+        "get_rotation_single"
+        help = "Get rotation of the ice floes"
+        action = :command
     end
 
     command_common_args = []

IFTPipeline.jl/src/h5.jl

@@ -139,7 +139,7 @@ function makeh5files(;
 end
 
 function makeh5files_single(;
-    passtime::DateTime,
+    passtime::ZonedDateTime,
     iftversion::Union{String,Nothing}=nothing,
     truecolor::String,
     falsecolor::String,
@@ -148,7 +148,7 @@ function makeh5files_single(;
     output::String,
 )
     latlondata = getlatlon(truecolor)
-    ptsunix = Int64(Dates.datetime2unix(passtime))
+    ptsunix = Int64(Dates.datetime2unix(DateTime(passtime)))
     labeled_ = load_labeled_img(labeled)
 
     if isnothing(iftversion)

IFTPipeline.jl/src/rotations.jl

@@ -0,0 +1,136 @@
+using DataFrames
+using TimeZones
+using Dates
+using CSV
+
+# Base.tryparse(::Type{ZonedDateTime}, str) = ZonedDateTime
+# default_format(::Type{ZonedDateTime}) = Format("yyyy-mm-dd\\THH:MM:SS.sZ")
+
+"""
+Make a CSV of pairwise rotations between floes detected on adjacent days. 
+
+Loads the floes from the `input` CSV file, and uses the columns:
+- `floe` ID
+- `satellite` name
+- `mask` – the binary mask (choose a column using argument `mask_column`)
+- `passtime` in ISO8601 format (with trailing Z or +00:00), e.g. 2022-09-11T09:21:00+00:00  (choose a column using argument `time_column`)
+- `date` of the overpass in YYYY-MM-DD format
+
+Returns a CSV with one row per floe comparison. 
+In the following, `i=1` means the earlier observation, `i=2` the later.
+
+Columns returned:
+- `ID` of the floe
+- Angle measures `theta_<deg,rad>` – angle between floe image in degrees or radians
+- Time measurements:
+  - `delta_time_sec` – number of seconds between overpass in the two measurements
+  - `omega_<deg,rad>_per_<sec,hour,day>` – mean angular velocity of rotation in degrees or radians per second hour or day.
+- Metadata 
+  - `satellite<i>` – which satellite measurement `i` was from
+  - `date<i>` – which date measurement `i` was taken
+  - `datetime<i>` – which UTC time measurement `i`'s overpass occurred
+- Original data
+  - `mask<i>` – the binary mask used for the measurement
+"""
+function get_rotation_single(;
+    input::String, output::String, mask_column=:mask, time_column=:passtime
+)
+    input_df = DataFrame(CSV.File(input))
+
+    input_df[!, mask_column] = eval.(Meta.parse.(input_df[:, mask_column]))
+    input_df[!, time_column] = ZonedDateTime.(input_df[:, time_column])
+
+    results = []
+    for row in eachrow(input_df)
+        append!( # adds the 0 – n measurements from `get_rotation_measurements` to the results array
+            results,
+            get_rotation_measurements(
+                row, input_df; mask_column=mask_column, time_column=time_column
+            ),
+        )
+    end
+    results_df = DataFrame(results)
+    @info results_df
+
+    FileIO.save(output, results_df)
+    return results_df
+end
+
+function get_rotation_measurements(
+    measurement::DataFrameRow, df::DataFrame; mask_column, time_column
+)
+    filtered_df = subset(
+        df,
+        :ID => ByRow(==(measurement[:ID])),
+        :date => ByRow(==(measurement[:date] - Dates.Day(1))),
+    )
+
+    results = [
+        get_rotation_measurements(
+            earlier_measurement, measurement; mask_column, time_column
+        ) for earlier_measurement in eachrow(filtered_df)
+    ]
+
+    return results
+end
+
+function get_rotation_measurements(
+    row1::DataFrameRow, row2::DataFrameRow; mask_column, time_column
+)
+    theta_rad = get_rotation(row1[mask_column], row2[mask_column])
+    theta_deg = rad2deg(theta_rad)
+
+    dt = row2[time_column] - row1[time_column]
+    dt_sec = dt / Dates.Second(1)
+    dt_hour = dt / Dates.Hour(1)
+    dt_day = dt / Dates.Day(1)
+
+    omega_deg_per_sec = (theta_deg) / (dt_sec)
+    omega_deg_per_hour = (theta_deg) / (dt_hour)
+    omega_deg_per_day = (theta_deg) / (dt_day)
+
+    omega_rad_per_sec = (theta_rad) / (dt_sec)
+    omega_rad_per_hour = (theta_rad) / (dt_hour)
+    omega_rad_per_day = (theta_rad) / (dt_day)
+
+    return (
+        ID=row1.ID,
+        theta_deg,
+        theta_rad,
+        delta_time_sec=dt_sec,
+        omega_deg_per_sec,
+        omega_deg_per_hour,
+        omega_deg_per_day,
+        omega_rad_per_sec,
+        omega_rad_per_hour,
+        omega_rad_per_day,
+        satellite1=row1.satellite,
+        satellite2=row2.satellite,
+        date1=row1.date,
+        date2=row2.date,
+        datetime1=row1[time_column],
+        datetime2=row2[time_column],
+        mask1=row1[mask_column],
+        mask2=row2[mask_column],
+    )
+end
+
+"""
+Get the angle in radians to rotate mask1 to mask2.
+"""
+function get_rotation(
+    mask1, mask2; mxshift::Tuple{Int64,Int64}=(100, 100), mxrot::Float64=Float64(pi)
+)
+    affine_map, _ = IceFloeTracker.Register.RegisterQD.qd_rigid(
+        IceFloeTracker.centered(mask1),
+        IceFloeTracker.centered(mask2),
+        mxshift,
+        mxrot;
+        print_interval=typemax(Int),
+    )
+    linear_map = affine_map.linear
+    cosθ = linear_map[1, 1]
+    sinθ = linear_map[2, 1]
+    θ = atan(sinθ, cosθ)
+    return θ
+end

IFTPipeline.jl/src/tracker.jl

@@ -42,7 +42,7 @@ end
     track_single(
         imgs::Array{String},
         props::Array{String},
-        passtimes::Array{DateTime},
+        passtimes::Array{ZonedDateTime},
         latlon::String,
         output::String,
         ...
@@ -65,7 +65,7 @@ Following are the default set of thresholds `condition_thresholds` used for floe
 function track_single(;
     imgs::Array{String},
     props::Array{String},
-    passtimes::Array{DateTime},
+    passtimes::Array{ZonedDateTime},
     latlon::String,
     output::String,
     dist::Array{Int}=[15, 30, 120],
@@ -139,7 +139,7 @@ function track_single(;
     adduuid!(props_)
 
     tracked_floes = long_tracker(props_, condition_thresholds, mc_thresholds)
-    FileIO.save(output, select!(tracked_floes, Not(:mask, :psi)))
+    FileIO.save(output, select!(tracked_floes, Not(:psi)))
     return tracked_floes
 end
 

IFTPipeline.jl/test/Project.toml

@@ -6,6 +6,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"