|
| 1 | +# # Fluxnet forcing data: LAI, radiation, and atmospheric variables |
| 2 | + |
| 3 | +# To access the forcing data (LAI from MODIS, SW_d, LW_d, T_air, q_air, |
| 4 | +# P_air, and precipitation from fluxtower data), you first need the |
| 5 | +# the fluxtower site ID. |
| 6 | +# Currently, ClimaLand provides an interface for working with four |
| 7 | +# fluxtower sites; adding support for a much larger set of sites is |
| 8 | +# in progress. The sites we support are Vaira Ranch (US_Var), |
| 9 | +# Missouri Ozark (US-MOz), Niwot Ridge (US-NR1), and Harvard Forest |
| 10 | +# (US-Ha1). |
| 11 | + |
| 12 | +using Dates |
| 13 | +import ClimaParams as CP |
| 14 | +using ClimaLand |
| 15 | +using ClimaLand.Domains: Column |
| 16 | +import ClimaLand.Parameters as LP |
| 17 | +using DelimitedFiles |
| 18 | +using CairoMakie |
| 19 | +import ClimaLand.FluxnetSimulations as FluxnetSimulations |
| 20 | + |
| 21 | +# Pick a site ID; convert the dash to an underscore: |
| 22 | +site_ID = "US-NR1"; |
| 23 | +site_ID_val = FluxnetSimulations.replace_hyphen(site_ID) |
| 24 | +# The functions we use below use multiple dispatch, treating the |
| 25 | +# site_ID_val as a Julia type. |
| 26 | + |
| 27 | +# First, we need the latitude and longitude of the site. These are used |
| 28 | +# to get the zenith angle as a function of time, and to look up |
| 29 | +# default parameters using the global ClimaLand parameter maps. We also |
| 30 | +# need the offset of the local time of the site in hours from UTC. This is |
| 31 | +# because ClimaLand simulations are carried out in UTC, while the fluxtower |
| 32 | +# data is reported in local time. |
| 33 | +(; time_offset, lat, long) = |
| 34 | + FluxnetSimulations.get_location(FT, Val(site_ID_var)); |
| 35 | + |
| 36 | +# ClimaLand also needs to know the height at which the atmospheric data |
| 37 | +# was recorded. |
| 38 | +(; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID)); |
| 39 | + |
| 40 | +# It is also useful to know the bounds of the data, |
| 41 | +# in UTC, to use as the start and stop date of the simulation. |
| 42 | +(start_date, stop_date) = |
| 43 | + FluxnetSimulations.get_data_dates(site_ID, time_offset); |
| 44 | + |
| 45 | +# Define the floating point precision desired (64 or 32 bit), and get the |
| 46 | +# parameter set holding constants used across CliMA Models. |
| 47 | +const FT = Float32; |
| 48 | +earth_param_set = LP.LandParameters(FT); |
| 49 | + |
| 50 | +# Now we can construct the forcing objects. Under the hood, this |
| 51 | +# function finds the local path to the fluxtower data (and downloads it |
| 52 | +# if it is not present) and reads the data in. It then creates |
| 53 | +# two objects, one called `atmos`, of type `PrescibedAtmosphere`, and |
| 54 | +# one called `radiation`, of type `PrescribedRadiativeFluxes`. These |
| 55 | +# encode the data in interpolating functions which allow us to |
| 56 | +# estimate the forcing at any time during the simulation using linear |
| 57 | +# interpolation across gaps in the data. |
| 58 | +(; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet( |
| 59 | + site_ID, |
| 60 | + lat, |
| 61 | + long, |
| 62 | + time_offset, |
| 63 | + atmos_h, |
| 64 | + start_date, |
| 65 | + earth_param_set, |
| 66 | + FT, |
| 67 | +); |
| 68 | +# The atmosphere object holds the air temperature, pressure, specific |
| 69 | +# humidity, wind speed, and liquid and solid precipitation fluxes. |
| 70 | +# Since many fluxtower sites do not measure the snow and rain fluxes |
| 71 | +# separately, we estimate them internally. This is optional, and by providing |
| 72 | +# the kwarg `split_precip = false`, you can change the behavior to |
| 73 | +# return all measured precip as a liquid water flux. |
| 74 | +# The radiation object holds the downwelling short and long wave |
| 75 | +# radiative fluxes. The diffuse fraction is estimated internally |
| 76 | +# using an empirical function, and the zenith angle is computed |
| 77 | +# analytically. |
| 78 | + |
| 79 | +# The simulation time is measured in seconds since the `start_date`. We |
| 80 | +# can get the atmospheric temperature at the `start_date` as follows: |
| 81 | +sim_time = 0.0 |
| 82 | +T = [NaN] |
| 83 | +evaluate!(T, atmos.T, sim_time); |
| 84 | +@show T |
| 85 | +# Note that `evaluate!` updates `T` in place. This is important: in our |
| 86 | +# simulations, we allocate memory for the forcing, and then update the values |
| 87 | +# at each step. If we did not do this, the dynamic memory allocation |
| 88 | +# would cause the simulation to be incredibly slow. |
| 89 | +# We can also plot the interpolated air temperature for the first day: |
| 90 | +sim_times = 0.0:1800.0:86400.0 # one day in seconds |
| 91 | +air_temps = [atmos.T(t) for t in sim_times] # allocates |
| 92 | +fig = CairoMakie.Figure() |
| 93 | +ax = CairoMakie.Axis( |
| 94 | + fig[1, 1], |
| 95 | + ylabel = "Temperature (K)", |
| 96 | + xlabel = "Date", |
| 97 | + title = "Near-surface air temperature at Niwot Ridge", |
| 98 | +) |
| 99 | +lines!(ax, sim_times .+ start_date, air_temps) |
| 100 | +CairoMakie.save("air_temp.png", fig); |
| 101 | +#  |
| 102 | + |
| 103 | +# We do something very similar with LAI, but now the data is coming |
| 104 | +# from MODIS. In this case, we linearly interpolate from the gridded |
| 105 | +# MODIS data to the latitude and longitude of the site. To do spatial |
| 106 | +# interpolation, we need to create a ClimaLand domain. |
| 107 | +domain = Column(; zlim = (-3.0, 0.0), nelements = 10, longlat = (long, lat)) |
| 108 | +surface_space = domain.space.surface |
| 109 | +# Get the paths to each year of MODIS data within the start and stop dates. |
| 110 | +modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths(; |
| 111 | + context = ClimaComms.context(surface_space), |
| 112 | + start_date, |
| 113 | + end_date = stop_date, |
| 114 | +) |
| 115 | +LAI = ClimaLand.prescribed_lai_modis( |
| 116 | + modis_lai_ncdata_path, |
| 117 | + surface_space, |
| 118 | + start_date, |
| 119 | +); |
| 120 | + |
| 121 | +# Just like with the air temperature, the LAI is an object that we can use |
| 122 | +# to linearly interpolate observed LAI to any simulation time. |
| 123 | + |
| 124 | +# It can also be useful to know the maximum LAI at a site. To do so, we |
| 125 | +# can call, for the first year of data (first element of `modis_lai_ncdata_path`): |
| 126 | +maxLAI = |
| 127 | + FluxnetSimulations.get_maxLAI_at_site(modis_lai_ncdata_path[1], lat, long) |
0 commit comments