Merge pull request #244 from JuliaClimate/hurst

new hurst code

Author: Balinus

Project.toml

@@ -6,21 +6,40 @@ version = "0.24.2"
 
 [deps]
 CFTime = "179af706-886a-5703-950a-314cd64e0468"
+CircularArrays = "7a955b69-7140-5f4e-a0ed-f168c5e2e749"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
+Extremes = "fe3fe864-1b39-11e9-20b8-1f96fa57382d"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LombScargle = "fc60dff9-86e7-5f2f-a8a0-edeadbb75bd9"
 LongMemory = "f5f8e4a8-cb56-40ea-a13e-090cc212d358"
 MarSwitching = "b2b8a7f1-da70-4e5f-beaa-e2774ec39c2f"
 NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c"
 Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99"
 
 [compat]
+CFTime = "0.1"
+CircularArrays = "1.3"
+DataFrames = "1.1"
+DimensionalData = "0.27, 0.28, 0.29"
+Extremes = "1"
+Interpolations = "0.13, 0.14, 0.15"
+LombScargle = "1"
+LongMemory = "0.1.2"
+MarSwitching = "0.2"
+NetCDF = "0.10, 0.11, 0.12"
+Polynomials = "4"
+Shapefile = "0.10, 0.11, 0.12, 0.13"
+StableRNGs = "1"
+Statistics = "1.11"
+YAXArrays = "0.6"
+Zarr = "0.9"
 julia = "1.10"
 
 [extras]

src/ClimateTools.jl

@@ -1,8 +1,10 @@
 module ClimateTools
 
+using CircularArrays
 using CFTime
 using Dates
 using DimensionalData
+using Extremes
 using LombScargle
 using LongMemory
 using Statistics
@@ -13,26 +15,37 @@ using Polynomials
 using Interpolations
 using DataFrames
 using MarSwitching
+using Shapefile
 
 include("aggregate.jl")
 include("autocorrelation.jl")
 include("biascorrect.jl")
 include("climatology.jl")
 include("ensembles.jl")
-#include("functions.jl")
+include("functions.jl")
+include("gev.jl")
 include("markov.jl")
 include("power.jl")
 include("plotting.jl")
 include("processERA5.jl")
+include("regrid.jl")
+include("spatial.jl")
+include("statistics.jl")
 include("utils.jl")
 
 
 export daily_fct, climato_tp, subsample, dates_builder_yearmonth, dates_builder_yearmonth_hardcode, dates_builder_yearmonthday, dates_builder_yearmonthday_hardcode, diff, cumsum, yearly_clim
+export ERA5Land_dailysum
+export m2mm
 export yearly_clim
+export quantiles
+export regrid_cube
 export qqmap, qqmap_bulk
 export ensemble_fct
 export autocorrelation
+export hurst
 export MSModel
+export rlevels_cube
 
 
 end

src/EFI.jl

@@ -0,0 +1,25 @@
+function estimate_prob(xout, xin1, xin2; da_cqt_quantile)
+    
+    # v_cqt = convert.(Float64,collect(da_cqt[31,31,:].data)) # the grid where we want to evaluate the interpolated value (ex. coords variable below)
+    # v_eqt = collect(da_eqt[:,31,31].data) # the x-coordinate of the initial data
+    
+    # v_cqt_quantile = collect(da_cqt.quantile); # The y-coordinates of the data points, same length as `xp`.
+    
+    itp = linear_interpolation(Interpolations.deduplicate_knots!(xin2), da_cqt_quantile, extrapolation_bc=Flat());
+    xout .= @. itp(xin1)
+    
+end
+
+
+function estimate_prob(da_cqt::YAXArray, da_eq::YAXArray)
+
+    Quantiles = collect(da_cqt.quantile)
+
+    # Dimensions
+    indims_da_cqt = InDims("quantile")
+    indims_da_eqt = InDims("Quantiles")
+    outdims = OutDims(Dim{:Probability}(Quantiles))
+
+    return mapCube(estimate_prob, (da_cqt, da_eq), indims=(indims_da_cqt, indims_da_eqt), outdims=outdims, da_cqt_quantile=Quantiles)
+    
+end

src/aggregate.jl

@@ -63,9 +63,9 @@ function daily_max_mean_min(cube::YAXArray, kwargs...)
     index_in_cube = [findall(==(i), time_index) for i in unique(time_index)]   
 
     # Dimensions
-    indims = InDims("Ti")    
+    indims = InDims("time")    
     # outdims = OutDims(RangeAxis("time", dates_builder_yearmonthday(new_dates)), CategoricalAxis("stats",["max","mean","min"]))
-    outdims = OutDims(Dim{:Ti}(dates_builder_yearmonthday(new_dates)), CategoricalAxis("stats",["max","mean","min"]))
+    outdims = OutDims(Dim{:time}(dates_builder_yearmonthday(new_dates)), CategoricalAxis("stats",["max","mean","min"]))
 
     mapCube(daily_max_mean_min, cube, indims=indims, outdims=outdims, index_list=index_in_cube)
 end
@@ -87,7 +87,7 @@ A new YAXArray cube with the computed percentiles along the "time", "longitude",
 function percentiles(cube::YAXArray, quantiles=[0.1, 0.5, 0.9];lonname="longitude", latname="latitude")    
 
     indims = InDims("number")    
-    outdims = OutDims("Ti", lonname, latname)
+    outdims = OutDims("time", lonname, latname)
 
     mapCube(percentiles, cube, indims=indims, outdims=outdims, quantiles=quantiles)
 end
@@ -124,8 +124,8 @@ Compute the difference between consecutive time steps in the given `cube`.
 function diff(cube::YAXArray)    
 
 
-    indims = InDims("Ti")    
-    outdims = OutDims("Ti")
+    indims = InDims("time")    
+    outdims = OutDims("time")
 
     mapCube(diff, cube, indims=indims, outdims=outdims)
 end
@@ -161,8 +161,8 @@ A new YAXArray cube with the cumulative sum computed along the "time" dimension.
 """
 function cumsum(cube::YAXArray, kwargs...)    
 
-    indims = InDims("Ti")    
-    outdims = OutDims("Ti")
+    indims = InDims("time")    
+    outdims = OutDims("time")
 
     mapCube(cumsum, cube, indims=indims, outdims=outdims)
 end
@@ -212,14 +212,14 @@ A new `YAXArray` with aggregated data on a daily basis.
 
 """
 function daily_fct(cube::YAXArray; fct::Function=mean, shifthour=0, kwargs...)
-    time_to_index = cube.Ti .+ Dates.Hour(shifthour)
+    time_to_index = cube.time .+ Dates.Hour(shifthour)
     time_index = yearmonthday.(time_to_index)
     new_dates = unique(time_index)
     index_in_cube = [findall(==(i), time_index) for i in unique(time_index)]   
 
     # Dimensions
-    indims = InDims("Ti")        
-    outdims = OutDims(Dim{:Ti}(dates_builder_yearmonthday(new_dates)))
+    indims = InDims("time")        
+    outdims = OutDims(Dim{:time}(dates_builder_yearmonthday(new_dates)))
 
     mapCube(daily_fct, cube, indims=indims, outdims=outdims, fct=fct, index_list=index_in_cube)
 end
@@ -254,14 +254,14 @@ A new YAXArray containing the yearly climatology.
 """
 
 function yearly_clim(cube::YAXArray; fct::Function=mean, kwargs...)
-    time_to_index = cube.Ti;
+    time_to_index = cube.time;
     time_index = year.(time_to_index);
     new_dates = unique(time_index);
     index_in_cube = [findall(==(i), time_index) for i in unique(time_index)]
 
     # Dimensions
-    indims = InDims("Ti")        
-    outdims = OutDims(Dim{:Ti}(dates_builder_yearmonthday_hardcode(new_dates, imois=7, iday=1)))
+    indims = InDims("time")        
+    outdims = OutDims(Dim{:time}(dates_builder_yearmonthday_hardcode(new_dates, imois=7, iday=1)))
 
     mapCube(yearly_clim, cube, fct=fct, indims=indims, outdims=outdims, index_list=index_in_cube)
 end
@@ -295,8 +295,8 @@ function daily_max(cube::YAXArray, kwargs...)
     index_in_cube = [findall(==(i), time_index) for i in unique(time_index)]   
 
     # Dimensions
-    indims = InDims("Ti")        
-    outdims = OutDims(Dim{:Ti}(dates_builder_yearmonthday(new_dates)))
+    indims = InDims("time")        
+    outdims = OutDims(Dim{:time}(dates_builder_yearmonthday(new_dates)))
 
     mapCube(daily_max, cube, indims=indims, outdims=outdims, index_list=index_in_cube)
 end

src/analysis.jl

@@ -19,9 +19,48 @@ The autocorrelation of the input dataset `ds` along the time dimension.
 function autocorrelation(ds; lags=30)
 
     # Dimensions
-    indims = InDims("Ti")
+    indims = InDims("time")
     outdims = OutDims(Dim{:lags}(collect(1:lags)))    
 
     return mapCube(ClimateTools.autocorrelation, ds, indims=indims, outdims=outdims, lags=lags, nthreads=Threads.nthreads())
 
+
+end
+
+function hurst(xout, xin; k::Int=20)
+    d = LongMemory.ClassicEstimators.rescaled_range_est(Array(xin), k=k)
+    xout .= d + 0.5
+    # xout .= LongMemory.autocorrelation(Array(xin), lags)   
+end
+
+
+
+"""
+    hurst(ds; k)
+
+Calculate the Hurst exponent for a time series of dataset `ds`. The Hurst exponent is a measure of the long-term memory of time series data. It can be used to classify time series as:
+
+- **0.5 < H < 1**: The time series is persistent, meaning that high values are likely to be followed by high values, and low values by low values.
+- **H = 0.5**: The time series is a random walk, indicating no correlation between values.
+- **0 < H < 0.5**: The time series is anti-persistent, meaning that high values are likely to be followed by low values, and vice versa.
+
+# Arguments
+- `ds`: YAXArray variable.
+
+# Returns
+- `H`: The Hurst exponent of the time series.
+
+
+"""
+function hurst(ds; k::Int=20)
+
+    # Dimensions
+    indims = InDims("time")
+    # outdims = OutDims(Dim{:hurst}(1))    
+    outdims = OutDims()    
+
+    return mapCube(ClimateTools.hurst, ds, indims=indims, outdims=outdims, k=k, nthreads=Threads.nthreads())
+
+
+
 end