Renaming some things

Author: Whyborn

src/domain/create_domain.jl

@@ -5,7 +5,6 @@ Domain{T, N, M}
 
 Data structure describing the computational domain for the simulation.
 
-$(FIELDS)
 """
 struct Domain{T, N, M}
     # The mask describing where the points are located on the grid.
@@ -27,56 +26,78 @@ struct Domain{T, N, M}
 end
 
 """
-Construct a Domain using vectors of longitudes and latitudes, a land_area_fraction array and a surface class mapping. The shape of the longitudes by latitudes should match the trailing dimensions of the land_area_fraction array. This signature is used for gridded domains.
+Construct a Domain using vectors of longitudes and latitudes, a vegetation_area_fraction array and a surface class mapping. The shape of the longitudes by latitudes should match the trailing dimensions of the vegetation_area_fraction array. This signature is used for gridded domains.
 
 If the lon_bnds and lat_bnds are not supplied, then they are inferred from the passed latitudes and longitudes. If `normalise` is true, the sum of fractions in a grid cell is normalised to 1.
 """
-function Domain(lons, lats, land_area_fraction{T, 3}, mapping; lon_bnds=nothing, lat_bnds=nothing, mask=nothing, normalise=false) where {T}
+function Domain(lons, lats, vegetation_area_fraction::Array{T, 3}, mapping; lon_bnds=nothing, lat_bnds=nothing, mask=nothing, normalise=false) where {T}
 
     # First validate the land area fraction. This includes checking the mapping is workable
     # with the passed land fractions i.e. the user hasn't denoted a tile index greater than
     # the number of slices in the array
-    @assert (size(land_area_fraction, 2), size(land_area_fraction, 3)) == (length(lons), length(lats)), "The size of the land_area_fraction and the supplied longitudes and latitudes are not compatible. Size of land_area_fraction: $(size(land_area_fraction)), with lons: $(length(lons)) and lats: $(length(lats))"
+    @assert size(vegetation_area_fraction, 1) == length(lons) && size(vegetation_area_fraction, 2) == length(lats) "The size of the vegetation_area_fraction and the supplied longitudes and latitudes are not compatible. Size of vegetation_area_fraction: $(size(vegetation_area_fraction)), with lons: $(length(lons)) and lats: $(length(lats))"
 
     # Convert the mapping values to vectors for consistent handling and array accessing,
     # then validate
     mapping = Dict(surface => inds isa Vector ? inds : [inds] for (surface, inds) in mapping)
-    validate_surface_mapping(land_area_fraction, mapping)
+    validate_surface_mapping(vegetation_area_fraction, mapping)
 
     # Were the bounds supplied?
     if (isnothing(lon_bnds))
         lon_bnds = infer_grid_bounds(lons)
     end
+    check_bounds_valid(lons, lon_bnds)
 
     if (isnothing(lat_bnds))
         lat_bnds = infer_grid_bounds(lats)
     end
+    check_bounds_valid(lats, lat_bnds)
+
+    # When creating the mask, we only want to work with the indices representing active surface
+    # types. So cut out a view of those indices from the complete land area fractions
+    active_inds = [ind for indices in values(mapping) for ind in indices]
+    active_area_fraction = @view(vegetation_area_fraction[:, :, active_inds])
+
+    # Now use this view to create the mask
+    if isnothing(mask)
+        # Note that non-land points may be treated as missing
+        frac_sums = dropdims(sum(active_area_fraction, dims=3), dims=3)
+        mask = @. !ismissing(frac_sums) && frac_sums > 0.0
+    else
+        @assert size(mask) == (size(vegetation_area_fraction, 1), size(vegetation_area_fraction, 2)), "Supplied mask does not match the shape of the land area fractions."
+    end
 
-    # Now we have handled everything that is specific to the gridded domains- reshape the
-    # data so that it can be treated in the same way as site runs.
-    # This means that the land_area_fraction should represent spatial points in a single
-    # dimension, and the mask should be a vector (if supplied)
-    land_area_fractions = reshape(land_area_fractions, (size(land_area_fractions), :))
-    mask = isnothing(mask) ? mask : vec(mask)
+    # Optionally normalise the fractions
+    if normalise
+        # Note that this modifies in place, in the variable that is a view into the
+        # original array
+        active_area_fraction .= active_area_fraction ./ sum(active_fractions, dims=2)
+    end
 
+    # Use the mask and land fractions to yield the index/fraction vectors
+    tiles, indices, fractions = process_vegetation_area_fractions(vegetation_area_fraction, mapping, mask)
+    
     # Vectorize the dimensions
     lons = repeat(lons, inner=length(lats), outer=1)
     lats = repeat(lats, inner=1, outer=length(lons))
     lon_bnds = repeat(lon_bnds, inner=(1, length(lats)), outer=(1, 1))
     lat_bnds = repeat(lat_bnds, inner=(1, 1), outer=(1, length(lons)))
 
     # Call the generic function used for both domain types
-    Domain(lons, lats, land_area_fractions, mapping; lon_bnds=lon_bnds, lat_bnds=lat_bnds, mask=mask, normalise=normalise)
+    Domain{eltype(vegetation_area_fraction), length(tiles), 2}(mask, lons, lats, lon_bnds, lat_bnds, tiles, indices, fractions)
 
 end
 
 """
-Construct a Domain using vectors of x and y coordinates (typically latitude and longitude), a land_area_fraction array and a surface class mapping. The longitudes and latitudes are paired, representing specific locations. This signature is used for site domains, and for processing of gridded domains after some preprocessing.
+Construct a Domain using vectors of x and y coordinates (typically latitude and longitude), a vegetation_area_fraction array and a surface class mapping. The longitudes and latitudes are paired, representing specific locations. This signature is used for site domains, and for processing of gridded domains after some preprocessing.
 
 If the lon_bnds and lat_bnds are not supplied, then it is assumed that the coordinates represent longitudes and latitudes, and the size of the grid cells is assumed to be 0.01 degrees.
 """
 
-function Domain(lons, lats, land_area_fraction{T 2}, mapping; lon_bnds=nothing, lat_bnds=nothing, mask=nothing, normalise=false) where {T}
+function Domain(lons, lats, vegetation_area_fraction::Matrix{T}, mapping; lon_bnds=nothing, lat_bnds=nothing, mask=nothing, normalise=false) where {T}
+
+    # Check that the lons and lats are valid
+    @assert length(lons) == length(lats) == size(vegetation_area_fraction, 2) "The longitudes, latitudes and number of columns in the vegetation_area_fraction array must be the same."
 
     # Create the lons and lats bounds if not supplied
     if isnothing(lon_bnds)
@@ -86,65 +107,70 @@ function Domain(lons, lats, land_area_fraction{T 2}, mapping; lon_bnds=nothing,
     if isnothing(lat_bnds)
         lat_bnds = infer_site_bounds(lats)
     end
-    # Verify that the bounds are valid if supplied
-    #
+
+    # Validate that the bounds are valid
+    check_bounds_valid(lons, lon_bnds)
+    check_bounds_valid(lats, lat_bnds)
+
     # Now we start processing the land fractions- we want to work with 2 copies of the
-    # underlying land_area_fractions. One which is the full array of data, that we can
+    # underlying vegetation_area_fractions. One which is the full array of data, that we can
     # index using the values in the mapping, and one which is a view of only the slices
     # included in the mapping.
-    # All the indices included in the mapping
     all_inds = [ind for indices in values(mapping) for ind in indices]
-    active_fractions = @view(land_area_fractions[all_inds, :])
+    active_area_fraction = @view(vegetation_area_fractions[:, all_inds])
 
     # Set the mask
     if (isnothing(mask))
         # Not supplied- created the mask. Take all of the "pages" of the land area array
         # that are mapped to a surface class, and check which locations have a non-zero
         # total surface fraction (and not missing which may represent ocean)
-        frac_sum = dropdims(sum(active_fractions, dims=1), dims=1)
+        frac_sum = dropdims(sum(active_area_fraction, dims=2), dims=2)
         mask = @. !ismissing(frac_sum) && frac_sum > 0.0
     else
         # Ensure that the mask size is compatible
-        @assert size(mask) == size(land_area_fraction, 2)
+        @assert length(mask) == size(vegetation_area_fraction, 1) "Supplied mask does not match the size of the vegetation_area_fraction."
     end
 
     # Optionally normalise the fractions
     if normalise
         # Note that this modifies in place, in the variable that is a view into the
         # original array
-        active_fractions .= active_fractions ./ sum(active_fractions, dims=1)
+        active_fractions .= active_area_fraction ./ sum(active_area_fraction, dims=2)
     end
 
     # Now create the vectors for the respective surface classes
-    tiles, indices, fractions = process_land_area_fractions_and_mapping(land_area_fraction, mapping, mask)
+    tiles, indices, fractions = process_vegetation_area_fractions_and_mapping(vegetation_area_fraction, mapping, vec(mask))
+
+    # Call the generic function used for both domain types
+    Domain{eltype(vegetation_area_fraction), length(tiles), 1}(mask, lons, lats, lon_bnds, lat_bnds, tiles, indices, fractions)
 
 end
 
 """
-Ensure that the specified mapping is valid with the given land area fractions. The only condition is that the mapping doesn't contain any indices outside the range available in the land area fractions i.e. 1 <= indx <= size(land_area_fraction, 1)
+Ensure that the specified mapping is valid with the given land area fractions. The only condition is that the mapping doesn't contain any indices outside the range available in the land area fractions i.e. 1 <= indx <= size(vegetation_area_fraction, 1)
 """
-function validate_surface_mapping(land_area_fraction, mapping)
+function validate_surface_mapping(vegetation_area_fraction, mapping)
     min_indx = 1
-    max_indx = size(land_area_fraction, 1)
+    max_indx = size(vegetation_area_fraction, 1)
     check_failed = false
     for (surface, inds) in mapping
         for ind in inds
             if !(min_indx <= ind <= max_indx)
-                @warn "Index $ind for surface class $surface falls outside the allowed range of 1 to size(land_area_fraction, 1), which is $max_indx"
+                @warn "Index $ind for surface class $surface falls outside the allowed range of 1 to size(vegetation_area_fraction, 1), which is $max_indx"
                 check_failed = true
             end
         end
     end
 
-    @assert !check_failed, "Supplied mapping is invalid for the given land area fractions" 
+    @assert !check_failed "Supplied mapping is invalid for the given land area fractions" 
 end
 
 """
 Infer the cell bounds for gridded domains, based on the list of cell centre locations passed.
 """
 function infer_grid_bounds(points::Vector{T}) where {T}
     # To infer the bounds, assume that the coordinates are ascending
-    @assert issorted(points), "Supplied points are not monotonically ascending"
+    @assert issorted(points) "Supplied points are not monotonically ascending"
 
     bounds = Matrix{T}(undef, 2, length(points))
     for i = 1:length(points)
@@ -158,8 +184,8 @@ function infer_grid_bounds(points::Vector{T}) where {T}
             lower = points[i] - points[i-1]
             upper = points[i+1] - points[i]
         end
-        bounds[1, i] = points[1] - 0.5 * lower
-        bounds[2, i] = points[1] + 0.5 * upper
+        bounds[1, i] = points[i] - 0.5 * lower
+        bounds[2, i] = points[i] + 0.5 * upper
     end
 
     bounds
@@ -168,7 +194,7 @@ end
 """
 Infer the cell bounds for site domains, based on the list of cell centre locations passed.
 """
-function infer_grid_bounds(points::Vector{T}) where {T}
+function infer_site_bounds(points::Vector{T}) where {T}
 
     bounds = Matrix{T}(undef, 2, length(points))
     bounds[1, :] = points .- 0.01
@@ -181,29 +207,32 @@ end
 Validate that the bounds are appropriate for the specified longitudes
 """
 function check_bounds_valid(points, point_bnds)
-    @assert all(@view(point_bnds[1, :]) .< points .< @view(point_bnds[2, :])), "Some bounds are do not bracket their associated point."
+    @assert all(@view(point_bnds[1, :]) .< points .< @view(point_bnds[2, :])) "Some bounds are do not bracket their associated point."
 end
 
 """
 Process the given land area fractions and the mapping to create the `tile`, `indices` and `fractions` components of the `Domain` derived type.
 """
-function process_land_area_fractions_and_mapping(land_area_fraction::Matrix{T}, mapping, mask::Vector) where {T}
+function process_vegetation_area_fractions(vegetation_area_fraction, mapping, mask)
 
+    # First ensure that the land area fraction and mask are vectorized
+    vegetation_area_fraction = reshape(vegetation_area_fraction, (:, size(vegetation_area_fraction, 1)))
+    mask = vec(mask)
 
     tiles = (); indices = (); fractions = ()
 
     for (surface_class, mapped_indices) in mapping
-        class_fractions = dropdims(sum(@view(land_area_fraction[mapped_indices]
-    components = Tuple(begin
-                           # We need to track the index in terms of land points
-                            land_counter=1
-                            for (ind, m) in enumerate(mask)
-                            if m
-
-                                
+        # Summate over the indices assigned to the surface class, then iterate through the non-zero and non-masked fractions
+        class_total = dropdims(sum(@view(vegetation_area_fraction[mapped_indices, :]), dims=2), dims=2)
+        indices_fracs = [(i, frac) for (i, (m, frac)) in enumerate(zip(mask, class_total)) if (m && frac > 0.0)]
 
-    for (surface_class, mapped_indices) in mapping
         tiles = (tiles..., surface_class)
+        indices = (indices..., getindex.(indices_fracs, 1))
+        fractions = (fractions..., getindex.(indices_fracs, 2))
+    end
+
+    tiles, indices, fractions
+end
 
 """
     define_on_domain(vars, domain::Domain{T}, dims=()
@@ -224,15 +253,3 @@ function define_on_domain(vars, domain::Domain{T}, dims=()) where {T}
 
     which_domain, ComponentVector(t)
 end
-
-function check_land_domain(domain::Domain, data::Array{T, 3}) where {T}
-    for i = 1:size(data, 3)
-        @assert all(domain.mask .== .!(ismissing.(@view(data[:, :, i])))), "Supplied data is expected to be on the same grid as the land mask, but it is not."
-    end
-end
-
-function check_land_domain(domain::Domain, data::Array{T, 4}) where {T}
-    for i = 1:size(data, 3), j = 1:size(data, 4)
-        @assert all(domain.mask .== .!(ismissing.(@view(data[:, :, i, j])))), "Supplied data is expected to be on the same grid as the land mask, but it is not."
-    end
-end

src/surface_types/ice_surface/define_ice_surface.jl

@@ -1,22 +1,22 @@
 # Not sure what traits a water surface may have yet
 
-const ICE_ALLOWED_TRAITS = Dict()
+const ICE_REQUIRED_TRAITS = Dict()
 
-# The master water surface type
-abstract type IceSurface <: SurfaceType end
+# The master ice surface type
+abstract type IceSurface <: SurfaceClass end
 
 """
     @IceType name, traits
 
 Define a new ice surface, which acts a functional type and subtype of WaterSurface.
 """
-macro IceType(name, ice_traits...)
-    traits = read_surface_traits(ICE_ALLOWED_TRAITS, ice_traits)
+macro IceSurface(name, ice_traits...)
+    traits = read_surface_traits(ICE_REQUIRED_TRAITS, ice_traits)
 
     esc(quote
         struct $(name) <: IceSurface end
 
         $name()
     end)
 end
-export @IceType, IceSurface
+export @IceSurface, IceSurface

src/surface_types/surface_types.jl

@@ -1,18 +1,21 @@
-abstract type SurfaceType end
+abstract type SurfaceClass end
 
-function read_surface_traits(allowed_traits, passed_traits)
+"""
+Read the surface traits passed, and check they form the complete set of traits required for the class.
+"""
+function read_surface_traits(required_traits, passed_traits)
     traits = Dict()
     for trait_spec in passed_traits
         if trait_spec.head == :(=)
-            if trait_spec.args[1] in keys(allowed_traits)
-                trait_options = allowed_traits[trait_spec.args[1]]
+            if trait_spec.args[1] in keys(required_traits)
+                trait_options = required_traits[trait_spec.args[1]]
                 if trait_spec.args[2] in keys(trait_options)
                     traits[trait_spec.args[1]] = trait_options[trait_spec.args[2]]
                 else
                     @error "$(trait_spec.args[2]) is an unrecognised value for the $(trait_spec.args[1]) trait. Options are $(keys(trait_options))."
                 end
             else
-                @error "$(trait_spec.args[1]) is an unrecognised PFT trait- must be one of $(keys(allowed_traits))"
+                @error "$(trait_spec.args[1]) is an unrecognised PFT trait- must be one of $(keys(required_traits))"
             end
         end
     end

src/surface_types/urban_surface/define_urban_surface.jl

@@ -1,22 +1,22 @@
 # Not sure what traits an urban surface may have
 
-const URBAN_ALLOWED_TRAITS = Dict()
+const URBAN_REQUIRED_TRAITS = Dict()
 
 # The master urban surface type
-abstract type UrbanSurface <: SurfaceType end
+abstract type UrbanSurface <: SurfaceClass end
 
 """
-    @UrbanType name, traits
+    @UrbanSurface name, traits
 
 Define a new Urban surface, which acts a functional type and subtype of WaterSurface.
 """
-macro UrbanType(name, urban_traits...)
-    traits = read_surface_traits(URBAN_ALLOWED_TRAITS, urban_traits)
+macro UrbanSurface(name, urban_traits...)
+    traits = read_surface_traits(URBAN_REQUIRED_TRAITS, urban_traits)
 
     esc(quote
         struct $(name) <: UrbanSurface end
 
         $name()
     end)
 end
-export @UrbanType, UrbanSurface
+export @UrbanSurface, UrbanSurface

src/surface_types/vegetated_surface/define_vegetated_surface.jl

@@ -9,20 +9,26 @@ const PHENOLOGY_TRAITS = Dict(:deciduous => Deciduous(),
                               )
 
 # The set of allowed traits
-const PFT_ALLOWED_TRAITS = Dict(:phenology => PHENOLOGY_TRAITS
+const VEGETATED_REQUIRED_TRAITS = Dict(:phenology => PHENOLOGY_TRAITS
                                 )
 
 # The master VegetationType
-abstract type VegetatedSurface <: SurfaceType end
+abstract type VegetatedSurface <: SurfaceClass end
 
 """
     @PFT name, traits
 
-Define a new PFT, which acts as a functional type and subtype of VegetationType, with
-the specified traits.
+An alias for @VegetatedSurface.
 """
 macro PFT(name, PFT_traits...)
-    traits = read_surface_traits(PFT_ALLOWED_TRAITS, PFT_traits)
+    :(@VegetatedSurface $name $PFT_traits)
+end
+
+"""
+Define a new VegetatedSurface, which acts as a functional type and subtype of VegetatedSurface, with the specified traits.
+"""
+macro VegetatedSurface(name, vegetated_traits...)
+    traits = read_surface_traits(VEGETATED_REQUIRED_TRAITS, vegetated_traits)
 
     esc(quote
         struct $(name) <: VegetatedSurface end
@@ -32,4 +38,4 @@ macro PFT(name, PFT_traits...)
         $name()
     end)
 end
-export @PFT, VegetatedSurface
+export @PFT, @VegetatedSurface, VegetatedSurface