From 836cd5920a6126a27b2477e1784d37e66198905d Mon Sep 17 00:00:00 2001 From: Joaquim Date: Tue, 13 May 2025 01:09:37 +0100 Subject: [PATCH] Fix NDWI and a couple other indices that were wrong. Needed to X -1 --- src/grid_at_sensor.jl | 4 +++- src/spectral_indices.jl | 7 +++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/grid_at_sensor.jl b/src/grid_at_sensor.jl index 9076679..320efca 100644 --- a/src/grid_at_sensor.jl +++ b/src/grid_at_sensor.jl @@ -167,6 +167,8 @@ function helper_find_sds(sds::String, info::String, ind_EOLs::Vector{UnitRange{I ((ind = findfirst(":lat", info)) === nothing) && (odeusse = true) elseif (startswith(sds, "lon")) ((ind = findfirst(":lon", info)) === nothing) && (odeusse = true) + else + odeusse = true end end odeusse && error("The band name -- " * sds * " -- does not exist") @@ -241,7 +243,7 @@ function get_lon_lat_qual(sds_lon::String, sds_lat::String, qual, inc) (inc[1] == 0.0) && (dx = diff(G.z[:, round(Int, size(G.z, 1)/2)]); dx = dx[dx .< 10]) lon = G.z[qual] G = gd2gmt(sds_lat); - (inc[2] == 0.0) && (dy = diff(G.z[round(Int, size(G.z, 2)/2), :]); dy = dy[dy .< 10]) + (inc[2] == 0.0) && (dy = diff(G.z[round(Int, size(G.z, 1)/2), :]); dy = dy[dy .< 10]) lat = G.z[qual] return lon, lat, dx, dy end diff --git a/src/spectral_indices.jl b/src/spectral_indices.jl index 34d892b..354c58e 100644 --- a/src/spectral_indices.jl +++ b/src/spectral_indices.jl @@ -693,15 +693,18 @@ function sp_indices(bnd1, bnd2, bnd3=nothing; index::String="", kwargs...) # NDBI, LSWI. Normalized difference water index. Gao 1996, Chen 2005; NDWI2 => (nir - swir2)/(nir + swir2) # Normalized difference red edge index. Gitelson and Merzlyak 1994; (redEdge2 - redEdge1)/(redEdge2 + redEdge1) # Normalized difference red edge index 2. Barnes et al 2000; (redEdge3 - redEdge1)/(redEdge3 + redEdge1) + + # Need to swap bnd1 and bnd2 for indices that are defined as (shorter_cdo - longer_cdo) + _b1, _b2 = (index == "MNDWI" || index == "NBRI" || index == "NDWI" || index == "NDWI2" || index == "NDWI2") ? (bnd2, bnd1) : (bnd1, bnd2) if (ismask) @inbounds Threads.@threads for k = 1:mn - t1 = bnd1[k]*i_tmax; t2 = bnd2[k]*i_tmax + t1 = _b1[k]*i_tmax; t2 = _b2[k]*i_tmax t = (t2 - t1) / (t1 + t2) (fcomp(t, threshold) && t <= 1) && (mask[k] = 255) end else @inbounds Threads.@threads for k = 1:mn - t1 = bnd1[k]*i_tmax; t2 = bnd2[k]*i_tmax + t1 = _b1[k]*i_tmax; t2 = _b2[k]*i_tmax t = (t2 - t1) / (t1 + t2) (t >= -1 && t <= 1) && (img[k] = t) end