Merge branch 'main' into albert-de-montserrat-patch-1

Author: boriskaus

.typos.toml

@@ -2,6 +2,7 @@
 MOR = "MOR"
 dum = "dum"
 Shepard = "Shepard"
+nin = "nin"
 
 [files]
 extend-exclude = ["tutorials/*.pvsm"]

Project.toml

@@ -1,14 +1,15 @@
 name = "GeophysicalModelGenerator"
 uuid = "3700c31b-fa53-48a6-808a-ef22d5a84742"
 authors = ["Boris Kaus", "Marcel Thielmann"]
-version = "0.7.4"
+version = "0.7.6"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+GDAL_jll = "a7073274-a066-55f0-b90d-d619367d196c"
 GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38"
 Geodesy = "0ef565a4-170c-5f04-8de2-149903a85f3d"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
@@ -26,6 +27,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+WhereTheWaterFlows = "ea906314-1493-4d22-a0af-f886a20c9fba"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [weakdeps]
@@ -43,7 +45,8 @@ Colors = "0.9 - 0.12"
 Downloads = "1"
 FFMPEG = "0.4"
 FileIO = "1"
-GLMakie = "0.10.3"
+GDAL_jll = "300.900.0 - 301.900.0"
+GLMakie = "0.10"
 GMT = "1"
 GeoParams = "0.2 - 0.6"
 Geodesy = "1"
@@ -60,6 +63,7 @@ Parameters = "0.9 - 0.12"
 SpecialFunctions = "1.0, 2"
 StaticArrays = "1"
 WriteVTK = "1"
+WhereTheWaterFlows = "0.10"
 julia = "1.9"
 
 [extras]
@@ -69,4 +73,4 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "GMT", "StableRNGs","GridapGmsh"]
+test = ["Test", "GMT", "StableRNGs", "GridapGmsh"]

ext/GLMakie_Visualisation.jl

@@ -3,7 +3,7 @@ module GLMakie_Visualisation
 
 using Statistics
 using GeophysicalModelGenerator: lonlatdepth_grid, GeoData, CartData, km, AbstractGeneralGrid
-import GeophysicalModelGenerator: visualise
+import GeophysicalModelGenerator: visualise, ustrip
 
 # We do not check `isdefined(Base, :get_extension)` as recommended since
 # Julia v1.9.0 does not load package extensions when their dependency is
@@ -14,7 +14,9 @@ else
     using ..GLMakie
 end
 
-export visualise
+import GLMakie: heatmap!, heatmap
+
+export visualise, heatmap, heatmap!
 
 println("Loading GLMakie extensions for GMG")
 
@@ -271,5 +273,121 @@ function visualise(Data::AbstractGeneralGrid; Topography=nothing, Topo_range=not
     return nothing
 end
 
+"""
+    heatmap(x::GeoData, args...; field=:Topography, kwargs...)
+heatmap for a 2D GeoData object (surface)
+"""
+function heatmap(x::GeoData, args...; field=:Topography, kwargs...)
+    @assert size(x.depth.val,3)==1
+    @assert hasfield(typeof(x.fields), field)
+
+    heatmap(x.lon.val[:,1], x.lat.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...)
+
+end
+
+"""
+    heatmap(x::GeoData, a::Array{_T,N}, args...; kwargs...)
+in-place heatmap for a 2D GeoData object (surface) with an array `a`. 
+"""
+function heatmap(x::GeoData, a::Array{_T,N}, args...; kwargs...) where{_T,N}
+    @assert size(x.depth.val,3)==1
+
+    if N==3
+        heatmap(x.lon.val[:,1], x.lat.val[1,:], ustrip.(a[:,:,1]), args...; kwargs...)
+    elseif N==2
+        heatmap(x.lon.val[:,1], x.lat.val[1,:], ustrip.(a), args...; kwargs...)
+    end
+
+end
+
+"""
+    heatmap(x::CartData, args...; field=:Topography, kwargs...)
+heatmap for a 2D CartData object (surface)
+"""
+function heatmap(x::CartData, args...; field=:Topography, kwargs...)
+    @assert size(x.z.val,3)==1
+    @assert hasfield(typeof(x.fields), field)
+
+    heatmap(x.x.val[:,1], x.y.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...)
+
+end
+
+"""
+    heatmap(x::CartData, a::Array{_T,N}, args...; kwargs...)
+in-place heatmap for a 2D CartData object (surface) with an array `a`. 
+"""
+function heatmap(x::CartData, a::Array{_T,N}, args...; kwargs...) where{_T,N}
+    @assert size(x.z.val,3)==1
+
+    if N==3
+        heatmap(x.x.val[:,1], x.y.val[1,:], ustrip.(a[:,:,1]), args...; kwargs...)
+    elseif N==2
+        heatmap(x.x.val[:,1], x.y.val[1,:], ustrip.(a), args...; kwargs...)
+    end
+
+end
+
+
+"""
+    heatmap!(x::GeoData, args...; field=:Topography, kwargs...)
+in-place heatmap for a 2D GeoData object (surface), 
+"""
+function heatmap!(x::GeoData, args...; field=:Topography, kwargs...)
+    @assert size(x.depth.val,3)==1
+    @assert hasfield(typeof(x.fields), field)
+
+    heatmap(x.lon.val[:,1], x.lat.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...)
+
+end
+
+"""
+    heatmap!(x::GeoData, a::Array{_T,N}, args...; kwargs...)
+in-place heatmap for a 2D GeoData object (surface) with an array `a`. 
+"""
+function heatmap!(x::GeoData, a::Array{_T,N}, args...; kwargs...) where{_T,N}
+    @assert size(x.depth.val,3)==1
+
+    if N==3
+        heatmap!(x.lon.val[:,1], x.lat.val[1,:], ustrip.(a[:,:,1]), args...; kwargs...)
+    elseif N==2
+        heatmap!(x.lon.val[:,1], x.lat.val[1,:], ustrip.(a), args...; kwargs...)
+    end
+
+end
+
+"""
+    heatmap!(x::CartData, args...; field=:Topography, colorbar=false, kwargs...)
+in-place heatmap for a 2D CartData object (surface)
+"""
+function heatmap!(x::CartData, args...; field=:Topography, colorbar=false, kwargs...)
+    @assert size(x.z.val,3)==1
+    @assert hasfield(typeof(x.fields), field)
+
+    data = ustrip.(x.fields[field][:,:,1])
+
+    fig,ax,hm = heatmap!(x.x.val[:,1], x.y.val[1,:], data, args...; kwargs...)
+
+    cb = nothing
+    if colorbar
+        cb = Colorbar(fig[1,2], limits=extrema(filter(!isnan,x.fields[field])), label=field)
+    end
+    
+    return fig_out
+end
+
+"""
+    heatmap!(x::CartData, a::Array{_T,N}, args...; kwargs...)
+in-place heatmap for a 2D CartData object (surface) with an array `a`. 
+"""
+function heatmap!(x::CartData, a::Array{_T,N}, args...; kwargs...) where{_T,N}
+    @assert size(x.z.val,3)==1
+
+    if N==3
+        heatmap!(x.x.val[:,1], x.y.val[1,:], ustrip.(a[:,:,1]), args...; kwargs...)
+    elseif N==2
+        heatmap!(x.x.val[:,1], x.y.val[1,:], ustrip.(a[:,:]), args...; kwargs...)
+    end
+
+end
 
 end

src/GeophysicalModelGenerator.jl

@@ -47,6 +47,8 @@ include("IO.jl")
 include("event_counts.jl")
 include("surface_functions.jl")
 include("movies_from_pics.jl")
+include("sea_lvl.jl")
+include("WaterFlow.jl")
 
 # Add optional routines (only activated when the packages are loaded)
 

src/WaterFlow.jl

@@ -0,0 +1,102 @@
+# These are routes to perform waterflow calculations (upstream area etc.) on a DEM
+
+using WhereTheWaterFlows
+import WhereTheWaterFlows: waterflows
+
+export waterflows
+
+"""
+    dlon, dlat = spacing(lon,lat)
+
+Computes the spacing with central differences
+"""
+function spacing(lon,lat)
+    dlon = zeros(size(lon.val)[1:2])
+    dlat = zeros(size(lat.val)[1:2])
+    @views dlon[2:end-1,:] = (lon.val[3:end,:,1] - lon.val[1:end-2,:,1])/2
+     dlon[1,:] = dlon[2,:]
+     dlon[end,:] = dlon[end-1,:]
+     dlat[:,2:end-1] = (lat.val[:,3:end,1] - lat.val[:,1:end-2,1])/2
+     dlat[:,1] = dlat[:,2]
+     dlat[:,end] = dlat[:,end-1]
+    
+    return dlon, dlat
+end
+
+"""
+    area_m2 = cell_area(Topo::GeoData)
+Returns the cell area for a Topographic dataset in m² (required for upstream area calculation)
+"""
+function cell_area(Topo::GeoData)
+
+    proj = ProjectionPoint(Lon=mean(Topo.lon.val[:]), Lat=mean(Topo.lat.val[:]))
+    Topo_cart = convert2CartData(Topo, proj)
+    dx, dy = spacing(Topo_cart.x, Topo_cart.y)
+
+    area_m2 = dx.*dy*1e6
+    return area_m2
+end
+
+
+"""
+    Topo_water, sinks, pits, bnds  = waterflows(Topo::GeoData; 
+        flowdir_fn=WhereTheWaterFlows.d8dir_feature, feedback_fn=nothing, drain_pits=true, bnd_as_sink=true,
+        rainfall = nothing,
+        minsize=300)
+
+Takes a GMG GeoData object of a topographic map and routes water through the grid. Optionally,
+you can specify `rainfall` in which case we accumulate the rain as specified in this 2D array instead of the cellarea. 
+This allows you to, for example, sum, up water if you have variable rainfall in the area.
+The other options are as in the `waterflows` function of the package `WhereTheWaterFlows`.
+
+Example
+===
+```julia
+# Download some topographic data
+julia> Topo = import_topo([6.5,7.3,50.2,50.6], file="@earth_relief_03s");
+
+# Flow the water through the area:
+julia> Topo_water, sinks, pits, bnds  = waterflows(Topo)
+julia> Topo_water
+GeoData 
+  size      : (961, 481, 1)
+  lon       ϵ [ 6.5 : 7.3]
+  lat       ϵ [ 50.2 : 50.59999999999999]
+  depth     ϵ [ 0.045 : 0.724]
+  fields    : (:Topography, :area, :slen, :dir, :nout, :nin, :c)
+
+```
+
+"""
+function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; feedback_fn=nothing, drain_pits=true, bnd_as_sink=true, rainfall=nothing, minsize=300) 
+
+    cellarea = cell_area(Topo)
+    cellarea_m2 = cellarea
+    if !isnothing(rainfall)
+        @assert typeof(rainfall) == Array{Float64,2}
+        cellarea = rainfall
+    end
+
+    dem = Topo.depth.val[:,:,1]
+
+    ni      = size(Topo.depth.val)
+    area    = zeros(ni)
+    slen    = zeros(Int64, ni)
+    dir     = zeros(Int8, ni)
+    nout    = zeros(Int8, ni)
+    nin     = zeros(Int8, ni)
+    c       = zeros(Int64, ni)
+
+    area[:,:,1], slen[:,:,1], dir[:,:,1], nout[:,:,1], nin[:,:,1], sinks, pits, c[:,:,1], bnds = waterflows(dem, cellarea, flowdir_fn;
+                        feedback_fn=feedback_fn, drain_pits=drain_pits, bnd_as_sink=bnd_as_sink)
+
+    catchment_large =  prune_catchments(c, minsize; val=0)
+    largest_catchment = catchment_large .== maximum(catchment_large)
+    largest_area = copy(area)
+    largest_area[.!largest_catchment] .= NaN
+    
+    log10_area = log10.(area)
+     
+    Topo_water =  addfield(Topo,(;area, slen, dir, nout, nin, c, cellarea_m2, catchment_large, log10_area, largest_catchment, largest_area))                        
+    return Topo_water, sinks, pits, bnds 
+end

src/data_types.jl

@@ -1464,14 +1464,6 @@ end
 extrema(d::FEData) = extrema(d.vertices, dims=2)
 size(d::FEData) = size(d.connectivity,2)
 
-"""
-    convert2FEData(d::Q1Data
-
-"""
-    convert2FEData(d::Q1Data)
-
-
-
 """
     X,Y,Z = coordinate_grids(Data::Q1Data; cell=false)
 
@@ -1531,4 +1523,4 @@ function convert2FEData(data::Q1Data)
     end
 
     return FEData(vertices,connectivity,  NamedTuple{keys(data.fields)}(data_fields),  NamedTuple{keys(data.cellfields)}(data_cellfields))
-end
+end

src/sea_level.jl

@@ -0,0 +1,48 @@
+export sea_level_files, SeaLevel, load_sea_level, curve_name
+
+const sea_level_path = "sea_level_data"
+
+const sea_level_files = Dict(
+    :Spratt_800ka => "Spratt2016-800ka.txt",
+    :Spratt_450ka => "Spratt2016-450ka.txt",
+    :Bintanja_3Ma => "Bintanja-3Ma.txt",
+    :Grant_153ka => "Grant2012_153ka.txt",
+    :Rohl_Bint_3Ma => "Rohl-Bint-3Ma.txt",
+    :Rohling2009_516ka => "Rohling2009-516ka.txt",
+    :Siddall2003_379ka => "Siddall2003-379ka.txt",
+    :Waelbroeck2002 => "Waelbroeck2002.txt",
+)
+
+struct SeaLevel{T}
+    elevation::Vector{T}
+    age::Vector{T}
+    name::Symbol
+
+    function SeaLevel(name::Symbol; flip_elevation = false, flip_age = false)
+        age, elevation = load_sea_level(
+            name; 
+            flip_age = flip_age,
+            flip_elevation = flip_elevation
+        )
+        new{eltype(age)}(age, elevation, name)
+    end
+end
+
+Base.getindex(x::SeaLevel, i::Int64) = x.elevation[i]
+Base.getindex(x::SeaLevel, i::Int64, j::Int64) = x.elevation[i], x.age[j]
+Base.getindex(x::SeaLevel, i::Tuple{Int64}) = x.elevation[i...], x.age[i...]
+Base.size(x::SeaLevel) = size(x.elevation)
+Base.eachindex(x::SeaLevel) = eachindex(x.elevation)
+Base.axes(x::SeaLevel) = axes(x.elevation)
+Base.length(x::SeaLevel) = length(x.elevation)
+curve_name(x::SeaLevel) = x.name
+
+function load_sea_level(name::Symbol; flip_elevation = false, flip_age = false)
+    fname = sea_level_files[name]
+    data = readdlm(joinpath(sea_level_path, fname))
+    h = data[:, 1]
+    age = data[:, 2]
+    flip_elevation && reverse!(h)
+    flip_age && reverse!(age)
+    return h, age
+end