enable gradients at subgrids too

src/parametrization/gradients.jl

@@ -22,24 +22,33 @@ lpars = localanisotropies(Gradients, grid, :CO2, 5)
 """
 function localanisotropies(::Type{Gradients}, obj, prop, window)
     # get dimensions
-    dims = size(domain(obj))
+    subd = domain(obj) isa SubDomain
+    dom = subd ? parent(domain(obj)) : domain(obj)
+    dims = size(dom)
     N = length(dims)
-    propv = Tables.getcolumn(Tables.columns(values(obj)), prop)
+    propv = rescale_min_max(Tables.getcolumn(Tables.columns(values(obj)), prop))
+
+    # extract gradients
+    img = if subd
+        im = fill(NaN, dims)
+        im[parentindices(domain(obj))] .= propv
+        im
+    else
+        reshape(propv, Size(dims))
+    end
 
-    # extract gradients (maybe need an extra method to deal with big datasets)
-    factor = maximum(propv) - minimum(propv)
-    factor = factor < 1000 ? 100 : 1
-    img = round.(Int, (factor .* reshape(propv, Size(dims)))) # temp solution
     g = imgradients(img, KernelFactors.sobel, "replicate")
 
     quat = Array{Quaternion}(undef, size(img)) # make some better way to store it
     m = Array{Vector}(undef, size(img)) # make some better way to store it
 
     Threads.@threads for i in CartesianIndices(img)
+        isnan(img[i]) && continue
         tensor = zeros(Float64, N, N)
         for ng in gridneighbors(img, i, window)
             for j in CartesianIndices((1:N, 1:N))
-                tensor[j] += g[j[1]][ng] * g[j[2]][ng]
+                val = g[j[1]][ng] * g[j[2]][ng]
+                !isnan(val) && (tensor[j] += val)
             end
         end
 
@@ -61,5 +70,24 @@ function localanisotropies(::Type{Gradients}, obj, prop, window)
         m[i] = λ / λ[1]
     end
 
-    LocalAnisotropy(vec(quat), reduce(hcat, vec(m)))
+    outid = findall(x -> !isnan(x), img)
+    LocalAnisotropy(vec(quat[outid]), reduce(hcat, vec(m[outid])))
+end
+
+
+function rescale_min_max(v)
+    # Extract non-NaN values
+    valid_vals = filter(!isnan, v)
+
+    # If there are no valid values, return the original vector
+    if isempty(valid_vals)
+        return v
+    end
+
+    # Compute min and max of non-NaN values
+    min_val = minimum(valid_vals)
+    max_val = maximum(valid_vals)
+
+    # Rescale the vector, keeping NaNs unchanged
+    return [isnan(x) ? NaN : (x - min_val) / (max_val - min_val) for x in v]
 end

src/parametrization/searchers.jl

@@ -4,7 +4,7 @@
 
 # grid neighborhood for gradients
 function gridneighbors(img, i::CartesianIndex, window::Int)
-    dim = Size(img)
+    dim = isa(img, StaticArray) ? Size(img) : size(img)
     minid(d) = max(1, i[d] - window)
     maxid(d) = min(dim[d], i[d] + window)
     idx = Tuple([minid(d):maxid(d) for d = 1:length(dim)])