Done data download prep for soil and dtm. TODO slope and aspect

forestmod · Jun 5, 2024 · b47669a · b47669a
1 parent 0bc0088
commit b47669a
Show file tree

Hide file tree

Showing 5 changed files with 359 additions and 384 deletions.
diff --git a/Project.toml b/Project.toml
@@ -4,7 +4,13 @@ authors = ["Antonello Lobianco <[email protected]>"]
 version = "0.1.0"
 
 [deps]
+ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
+BetaML = "024491cd-cc6b-443e-8034-08ea7eb7db2b"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+FTPClient = "01fcc997-4f28-56b8-8a06-30002c134abb"
 GenFSM_resources = "731fd045-8cb1-4036-8c04-a2a477f3265a"
-
-#[sources]
-#GenFSM_resources = {url = "[email protected]:forestmod/GenFSM_resources.jl.git", rev = "main"}
+Geomorphometry = "714e3e49-7933-471a-9334-4a6a65a92f36"
+RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
+Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
+ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
diff --git a/src/GenFSM_resources_init_fr.jl b/src/GenFSM_resources_init_fr.jl
@@ -6,7 +6,15 @@ Initializer for the France side of the resource module
 """
 module GenFSM_resources_init_fr
 
-export data_path, init!
+import Downloads 
+import ArchGDAL, Rasters, ZipFile, DataStructures, FTPClient
+import Geomorphometry # for slope and aspect 
+
+#using FTPClient, ZipFile, DataStructures #, RasterDataSources
+
+export data_path, init!!
+
+include("Utils.jl")
 
 #using Statistics, Random, Downloads # Standard Library (shipped with Julia)
 #using FTPClient
@@ -24,11 +32,12 @@ export data_path, init!
 
 #data_path   = joinpath(@__DIR__,"data")
 #ENV["RASTERDATASOURCES_PATH"] = data_path
-#include("getdata.jl")
 
-function init!(pixels,mask,settings)
-    println("hello w")
-    #temp_ rel_temp_output
+#ENV["RASTERDATASOURCES_PATH"] = "/home/lobianco/CloudFiles/beta-lorraine-sync/MrFOR/dev_GenFSM/cache/res/fr"
+
+include("Get_data.jl")
+
+function init!!(pixels,settings,mask)
     temp_path = joinpath(settings["temp_path"],"res","fr")
     cache_path = joinpath(settings["cache_path"],"res","fr")
     output_path = joinpath(settings["output_path"],"res","fr")
@@ -38,6 +47,10 @@ function init!(pixels,mask,settings)
     isdir(temp_path) || mkpath(temp_path)
     isdir(cache_path) || mkpath(temp_path)
     isdir(output_path) || mkpath(output_path)
+
+    get_data!(settings,mask)
+    println(settings)
+
     # Download the data:
     #- administrative for the region
     #- soil 
@@ -48,4 +61,9 @@ end
 
 
 
+
+
+
+
+
 end # module GenFSM_resources_init_fr
diff --git a/src/Get_data.jl b/src/Get_data.jl
@@ -0,0 +1,186 @@
+# See also Copernicus or CHELSA for climate:
+# https://cds.climate.copernicus.eu/cdsapp#!/search?type=dataset
+# https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=overview
+# https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/EUR11/documentation/CHELSA_EUR11_technical_documentation.pdf
+
+
+function get_data!(settings,mask)
+    input_rasters = Dict{Any,Any}()
+    input_rasters = merge(input_rasters,get_dtm(settings,mask))
+    input_rasters = merge(input_rasters,get_soil_data(settings,mask))
+    settings["res"]["fr"]["input_rasters"] = input_rasters 
+end
+
+# ------------------------------------------------------------------------------
+# Specific download functions
+
+function get_dtm(settings,mask)
+    data_path = joinpath(settings["res"]["fr"]["cache_path"],"dtm")
+    isdir(data_path) || mkpath(data_path)
+    force     = "dtm" in settings["res"]["fr"]["force_download"]
+    url       = settings["res"]["fr"]["data_sources"]["dtm_url"]
+    to        = mask
+    filename  = basename(url)
+    verbosity = settings["verbosity"]
+    verbose   = verbosity in ["HIGH", "FULL"] 
+    zip_destpath = joinpath(data_path,filename)
+    tif_destpath    = replace(zip_destpath,".zip" => ".tif")
+    tif_destpath_reg = replace(tif_destpath,".tif" => "_reg.tif")
+    dtm_var   = DataStructures.OrderedDict("dtm"=>tif_destpath_reg)
+    if ( isfile(tif_destpath) && (!force))
+        return dtm_var
+    end
+    verbose && @info "Downloading dtm file..."
+    Downloads.download(url,zip_destpath, verbose=verbose)
+    unzip(zip_destpath,data_path)
+    dtm_w   = Rasters.Raster(tif_destpath)
+    dtm_w   = Rasters.crop(rnge; to=to)
+
+    dtm_reg = Rasters.resample(dtm_w,to=to,method=:average)
+    write(tif_destpath_reg, dtm_reg,force=true)
+    rm(zip_destpath)
+    #rm(tif_destpath_reg)
+    return dtm_var
+end
+
+
+
+function get_soil_data(settings,mask)
+    #=
+    All soil datasets:
+    https://esdac.jrc.ec.europa.eu/resource-type/european-soil-database-soil-properties
+
+    Soil_physical 
+    https://esdac.jrc.ec.europa.eu/content/topsoil-physical-properties-europe-based-lucas-topsoil-data
+    https://www.sciencedirect.com/science/article/pii/S0016706115300173
+
+    Soil chemistry
+    https://esdac.jrc.ec.europa.eu/content/chemical-properties-european-scale-based-lucas-topsoil-data
+    https://www.sciencedirect.com/science/article/pii/S0016706119304768
+
+    Soil other derived data
+    https://esdac.jrc.ec.europa.eu/content/european-soil-database-derived-data
+    Hiederer, R. 2013. Mapping Soil Properties for Europe - Spatial Representation of Soil Database Attributes. Luxembourg: Publications Office of the European Union - 2013 - 47pp. EUR26082EN Scientific and Technical Research series, ISSN 1831-9424, doi:10.2788/94128
+    Hiederer, R. 2013. Mapping Soil Typologies - Spatial Decision Support Applied to European Soil Database. Luxembourg: Publications Office of the European Union - 2013 - 147pp. EUR25932EN Scientific and Technical Research series, ISSN 1831-9424, doi:10.2788/8728
+    =#
+
+    soil_path      = joinpath(settings["res"]["fr"]["cache_path"],"soil")
+    isdir(soil_path) || mkpath(soil_path)
+    soil_ph_url    = settings["res"]["fr"]["data_sources"]["soil_ph_url"]
+    soil_chem_url  = settings["res"]["fr"]["data_sources"]["soil_chem_url"]
+    soil_oth_url   = settings["res"]["fr"]["data_sources"]["soil_oth_url"]
+    soil_ph_vars   = settings["res"]["fr"]["data_sources"]["soil_ph_vars"]
+    soil_chem_vars = settings["res"]["fr"]["data_sources"]["soil_chem_vars"]
+    soil_oth_vars  = settings["res"]["fr"]["data_sources"]["soil_oth_vars"]
+    force          = "soil" in settings["res"]["fr"]["force_download"]
+    verbosity      = settings["verbosity"]
+    verbose        = verbosity in ["HIGH", "FULL"] 
+    to             = mask
+    soil_vars      = DataStructures.OrderedDict{String,String}()
+    n_soil_ph_vars = length(soil_ph_vars)
+    soil_texture_n_classes = settings["res"]["fr"]["data_sources"]["soil_texture_n_classes"]
+    crs = convert(Rasters.WellKnownText, Rasters.EPSG(settings["simulation_region"]["cres_epsg_id"]))
+
+    # Soil and chemistry variables....
+    for (i,var) in enumerate(vcat(soil_ph_vars,soil_chem_vars))
+        if i <= n_soil_ph_vars
+            urlname   = replace(soil_ph_url,"\${VAR}" => var)
+        else
+            urlname   = replace(soil_chem_url,"\${VAR}" => var)
+        end
+        zipname   = joinpath(soil_path,"$(var).zip")
+        zipfolder = joinpath(soil_path,var)
+        final_file = joinpath(soil_path,"$(var)_reg.tif") 
+        if (ispath(final_file) && !force )
+            soil_vars[var] = final_file
+            continue
+        end
+        Downloads.download(urlname,zipname, verbose=verbose)
+        unzip(zipname,zipfolder)
+        dir_files = readdir(zipfolder)
+        tif_file = dir_files[findfirst(x->match(r".*\.tif$", x) !== nothing, dir_files)]
+        saved_file = joinpath(zipfolder,tif_file)
+        orig_raster = Rasters.Raster(saved_file)
+        resampled_raster = Rasters.resample(orig_raster,to=to,method=:average)
+        write(final_file, resampled_raster, force=true)
+        soil_vars[var] = final_file
+        rm(zipname)
+        rm(zipfolder,recursive=true)
+    end
+
+    # Special: for TextureUSDA we extract the class as boolean values, as this is categorical
+    soil_texture_classes = 1:soil_texture_n_classes
+    texture_filename = joinpath(soil_path,"TextureUSDA_reg.tif")
+    texture_classes = expand_classes(texture_filename,soil_texture_classes;verbose=verbose,force=force,to=nothing)
+    delete!(soil_vars,"TextureUSDA")
+    soil_vars = DataStructures.OrderedDict(soil_vars..., texture_classes...)
+
+
+    # Other soil variables
+    urlname   = soil_oth_url
+    zipname   = joinpath(soil_path,basename(urlname))
+    zipfolder = joinpath(soil_path, split(basename(urlname),".")[1])
+    finalfolder = joinpath(soil_path,"soil_oth_vars")
+    if (ispath(finalfolder) && (!force) )
+        return soil_vars
+    end
+    ispath(finalfolder) || mkpath(finalfolder)
+    Downloads.download(urlname,zipname, verbose=verbose)
+    unzip(zipname,zipfolder)
+    for var in soil_oth_vars
+        saved_file = joinpath(zipfolder,"$(var).rst")
+        final_file = joinpath(finalfolder,"$(var)_reg.tif")
+        orig_raster = Rasters.Raster(saved_file,crs=crs) # ,missingval=0.0f0)  missing value is zero, but zero is used also as true value! 
+        resampled_raster = Rasters.resample(orig_raster,to=to,method=:average)
+        write(final_file, resampled_raster,force=true)
+        soil_vars[var] = final_file
+    end
+    rm(zipname)
+    rm(zipfolder,recursive=true)
+    return soil_vars
+end
+
+
+function get_ecoregions_data(data_path;force=false,verbose=false,to=nothing)
+    # landing page: https://www.eea.europa.eu/en/datahub/datahubitem-view/c8c4144a-8c0e-4686-9422-80dbf86bc0cb?activeAccordion=
+    #force   = false
+    #verbose = true
+    #to      = nothing
+    #to      = Yraster
+
+    subpath     = joinpath(data_path,"ecoregions")
+    classes     = collect([1:9;11:13]) # classes 10, 14 and 15 are out of the study area
+    urlname     = "https://sdi.eea.europa.eu/datashare/s/sTnNeQK69iYNgCe/download"
+    dataname    = "eea_r_3035_1_km_env-zones_p_2018_v01_r00"
+    zippath     = joinpath(subpath,"$(dataname).zip")
+    zipfolder   = joinpath(subpath,dataname)
+    tifname     = joinpath(zipfolder,"$(dataname).tif")
+    tifname_1km = joinpath(zipfolder,"ecoregions_1km.tif")
+    tifname_8km = joinpath(zipfolder,"ecoregions_8km.tif")
+    (ispath(zipfolder) && !force ) && return Dict("ecoregions" => tifname_8km)  #expand_classes(tifname_1km,classes;verbosity=verbose,force=false,to=to)
+
+    Downloads.download(urlname,zippath, verbose=verbose)
+    rm(zipfolder, force=true,recursive=true)
+    unzip(zippath,subpath) # the zip archive contain itself the folder with the data inside
+    rm(zippath)
+
+    # before the extraction of the classes we get the base raster res from 100m to 1km 
+    base_raster = Raster(tifname)
+    base_raster_1km = resample(base_raster, res=1000)
+    base_raster_8km = resample(base_raster_1km, to=to, method=:mode)
+    write(tifname_8km, base_raster_8km, force=true)
+    return Dict("ecoregions" => tifname_8km)
+
+    #=
+    #base_raster_1km_Float32 = Float32.(base_raster_1km)
+    # replacing 0 as missing value to 99, as we need 0 as a value !
+    missval       = missingval(base_raster_1km) 
+    base_raster_1km = map(x-> ( (x == missval) ? UInt8(99) : x), base_raster_1km)
+    base_raster_1km = rebuild(base_raster_1km; missingval=UInt8(99))
+    write(tifname_1km, base_raster_1km, force=true)
+
+    ecoregions_classes = expand_classes(tifname_1km,classes;verbosity=2,force=force,to=to)
+
+    return ecoregions_classes
+    =#
+end