From fd66967c877fd2aca521ddda55c749e502035b7f Mon Sep 17 00:00:00 2001 From: Dongdong Kong Date: Sat, 12 Oct 2024 14:32:38 +0800 Subject: [PATCH] third --- .gitignore | 2 + data/ee_extract.py | 74 ++++++++++++++++++++ data/ee_extractor.py | 54 +++----------- data/s1_GEE.qmd | 2 +- src/DataType/Constant.jl | 2 +- src/photosynthesis.jl | 147 ++++++++++++++++++++++++--------------- 6 files changed, 177 insertions(+), 104 deletions(-) create mode 100644 data/ee_extract.py diff --git a/.gitignore b/.gitignore index bb67a3e..71faa73 100644 --- a/.gitignore +++ b/.gitignore @@ -27,3 +27,5 @@ docs/reference .ipynb_checkpoints __pycache__ +data +examples diff --git a/data/ee_extract.py b/data/ee_extract.py new file mode 100644 index 0000000..e4b464f --- /dev/null +++ b/data/ee_extract.py @@ -0,0 +1,74 @@ +import argparse +import ee +import pandas as pd + +# Dongdong Kong +# js版本已经测试通过 +# 这是GPT3.5自动翻译的代码，可能会存在bug。 +# https://code.earthengine.google.com/027ef28c2f7f304be6c7ebed95f5e1cc?noload=1 + +# Initialize the Earth Engine library +ee.Initialize() + + +def image_extract_points(img, points, proj=None, set_date=False): + if proj == None: + proj = img.projection().getInfo() + options = { + 'collection': points, + 'reducer': 'first', + 'crs': proj['crs'], + 'crsTransform': proj['transform'], + 'tileScale': 16 + } + date = ee.Date(img.get('system:time_start')).format('yyyy-MM-dd') + + def tidy_props(f): + f = ee.Feature(None).copyProperties(f) + if set_date: + f = f.set('date', date) + return f + return img.reduceRegions(**options).map(tidy_props) + + +def col_extract_points(col, task, points): + proj = col.first().select(0).projection().getInfo() + res = col.map(lambda img: image_extract_points( + img, points, proj, set_date=True)).flatten() + + ee.batch.Export.table.toDrive( + collection=res, + description=task, + folder="gee", + fileNamePrefix=task, + fileFormat="GeoJSON" + ).start() + + +def export_image(img, task, prefix=""): + task = prefix + task + proj = img.projection().getInfo() + + # Export the image to Drive as a GeoTIFF file + ee.batch.Export.image.toDrive( + image=img, + description=task, + folder="gee", + crs=proj['crs'], + crsTransform=proj['transform'], + maxPixels=1e13, + fileNamePrefix=task, + fileFormat='GeoTIFF' + ).start() + + +def shp2df(fc, outfile=None): + data = fc.getInfo().get('features') + props = [x["properties"] for x in data] + df = pd.DataFrame(props) + if None != outfile: + df.to_csv(outfile, index=False) + return df + +def hello(): + print("Hello, World!") diff --git a/data/ee_extractor.py b/data/ee_extractor.py index 98560cf..b7cfdb8 100644 --- a/data/ee_extractor.py +++ b/data/ee_extractor.py @@ -1,3 +1,12 @@ +# %% +## add pwd to path +# import os +# import sys +# sys.path.append(os.getcwd()) + +from ee_extract import * + + # %% import argparse import ee @@ -13,23 +22,6 @@ # %% # Define the function to extract points from an image (internal) -def image_extract_points(img, points, proj=None, set_date=False): - if proj == None: proj = img.projection().getInfo() - options = { - 'collection': points, - 'reducer': 'first', - 'crs': proj['crs'], - 'crsTransform': proj['transform'], - 'tileScale': 16 - } - date = ee.Date(img.get('system:time_start')).format('yyyy-MM-dd') - - def tidy_props(f): - f = ee.Feature(None).copyProperties(f) - if set_date: f = f.set('date', date) - return f - - return img.reduceRegions(**options).map(tidy_props) def shp2df(fc, outfile=None): data = fc.getInfo().get('features') @@ -57,34 +49,6 @@ def shp2df(fc, outfile=None): # %% # Define the function to export an image to Drive -def export_image(img, task, prefix=""): - task = prefix + task - proj = img.projection().getInfo() - - # Export the image to Drive as a GeoTIFF file - ee.batch.Export.image.toDrive( - image=img, - description=task, - folder="gee", - crs=proj['crs'], - crsTransform=proj['transform'], - maxPixels=1e13, - fileNamePrefix=task, - fileFormat='GeoTIFF' - ).start() - -# Define the function to extract points from an image collection -def col_extract_points(col, task, points): - proj = col.first().select(0).projection().getInfo() - res = col.map(lambda img: image_extract_points(img, points, proj, set_date=True)).flatten() - - ee.batch.Export.table.toDrive( - collection=res, - description=task, - folder="gee", - fileNamePrefix=task, - fileFormat="GeoJSON" - ).start() def main(): parser = argparse.ArgumentParser(description="Google Earth Engine Command Line Tool") diff --git a/data/s1_GEE.qmd b/data/s1_GEE.qmd index dc0b0c8..9bf3619 100644 --- a/data/s1_GEE.qmd +++ b/data/s1_GEE.qmd @@ -3,7 +3,7 @@ ```{julia} using EarthEngine -# Initialize() +Initialize() dem = EE.Image("USGS/SRTMGL1_003") xy = Point(86.9250, 27.9881) diff --git a/src/DataType/Constant.jl b/src/DataType/Constant.jl index 265b4e8..52ed775 100644 --- a/src/DataType/Constant.jl +++ b/src/DataType/Constant.jl @@ -94,7 +94,7 @@ const ekc::Float64 = 80500.0 # Activation energy for K of CO2; J mol-1 const eko::Float64 = 14500.0 # Activation energy for K of O2, J mol-1 const erd::Float64 = 38000.0 # activation energy for dark respiration, eg Q10=2 const ektau::Float64 = -29000.0 # J mol-1 (Jordan and Ogren, 1984) -const tk_25::Float64 = 298.16 # absolute temperature at 25 C +const TK25::Float64 = 298.16 # absolute temperature at 25 C const toptvc::Float64 = 301.0 # optimum temperature for maximum carboxylation const toptjm::Float64 = 301.0 # optimum temperature for maximum electron transport const eabole::Float64 = 45162.0 # activation energy for bole respiration for Q10 = 2.02 diff --git a/src/photosynthesis.jl b/src/photosynthesis.jl index 862e717..6ea7461 100644 --- a/src/photosynthesis.jl +++ b/src/photosynthesis.jl @@ -1,5 +1,6 @@ function photosynthesis(Tc_old::Leaf, R::Leaf, Ci_old::Leaf, leleaf::Leaf, - Ta::Cdouble, ea::Cdouble, f_soilwater::Cdouble, b_h2o::Cdouble, m_h2o::Cdouble, + Ta::Cdouble, ea::Cdouble, f_soilwater::Cdouble, + g0_h2o::Cdouble, g1_h2o::Cdouble, Gb_o::Cdouble, Gb_u::Cdouble, Vcmax_sunlit::Cdouble, Vcmax_shaded::Cdouble, # output Gs_new::Leaf, Ac::Leaf, Ci_new::Leaf; version="c") @@ -11,30 +12,58 @@ function photosynthesis(Tc_old::Leaf, R::Leaf, Ci_old::Leaf, leleaf::Leaf, end Gs_new.o_sunlit, Ac.o_sunlit, Ci_new.o_sunlit = - fun(Tc_old.o_sunlit, R.o_sunlit, ea, Gb_o, Vcmax_sunlit, f_soilwater, b_h2o, m_h2o, + fun(Tc_old.o_sunlit, R.o_sunlit, ea, Gb_o, Vcmax_sunlit, f_soilwater, g0_h2o, g1_h2o, Ci_old.o_sunlit, Ta, leleaf.o_sunlit) Gs_new.o_shaded, Ac.o_shaded, Ci_new.o_shaded = - fun(Tc_old.o_shaded, R.o_shaded, ea, Gb_o, Vcmax_shaded, f_soilwater, b_h2o, m_h2o, + fun(Tc_old.o_shaded, R.o_shaded, ea, Gb_o, Vcmax_shaded, f_soilwater, g0_h2o, g1_h2o, Ci_old.o_shaded, Ta, leleaf.o_shaded) Gs_new.u_sunlit, Ac.u_sunlit, Ci_new.u_sunlit = - fun(Tc_old.u_sunlit, R.u_sunlit, ea, Gb_u, Vcmax_sunlit, f_soilwater, b_h2o, m_h2o, + fun(Tc_old.u_sunlit, R.u_sunlit, ea, Gb_u, Vcmax_sunlit, f_soilwater, g0_h2o, g1_h2o, Ci_old.u_sunlit, Ta, leleaf.u_sunlit) Gs_new.u_shaded, Ac.u_shaded, Ci_new.u_shaded = - fun(Tc_old.u_shaded, R.u_shaded, ea, Gb_u, Vcmax_shaded, f_soilwater, b_h2o, m_h2o, + fun(Tc_old.u_shaded, R.u_shaded, ea, Gb_u, Vcmax_shaded, f_soilwater, g0_h2o, g1_h2o, Ci_old.u_shaded, Ta, leleaf.u_shaded) end +""" +Ball-Berry stomatal conductance model + +# Arguments +- `A` : net photosynthesis rate (umol m-2 s-1) +- `Cs` : CO2 concentration at leaf surface (ppm) +- `RH` : relative humidity at leaf surface (0-1) +- `g0` : the minimum stomatal conductance (umol m-2 s-1) + +# Returns +- `gs` : stomatal conductance of co2 (umol m-2 s-1) +""" +function stomatal_conductance(A, Cs, RH, g0, g1, β_soil=1.0) + gs = g0 + β_soil * g1 * A * RH / Cs + return gs +end + # 这里依赖的变量过多，不易核对 """ +# Arguments +- `ea`: [kPa] +- `gb_w`: leaf laminar boundary layer conductance to H2O, [s m-1] +- `Vcmax25`: the maximum rate of carboxylation of Rubisco at 25 deg C, [umol m-2 s-1] + +- `cii` : intercellular CO2 concentration (ppm) +- `g0_h2o`: the minimum stomatal conductance to H2O, [mol m-2 s-1] +- `g1_h2o`: the slope of the stomatal conductance to H2O, [mol m-2 s-1] + +- `LH_leaf`: latent heat of vaporization of water at the leaf temperature, [W m-2] + # Intermediate variables -- `g_lb_c` : leaf laminar boundary layer condunctance to CO2 (mol m-2 s-1) +- `gb_c` : leaf laminar boundary layer condunctance to CO2 (mol m-2 s-1) - `RH_leaf` : relative humidity at leaf surface (0-1) - `gs_co2_mole`: stomatal conductance to CO2 (mol m-2 s-1) - `gs_h2o_mole`: stomatal conductance to h2o (mol m-2 s-1) @@ -50,69 +79,72 @@ end - `tau` : the specifity of Rubisco for CO2 compared with O2 - `Jₓ` : the flux of electrons through the thylakoid membrane (umol m-2 s-1) """ -@fastmath function photosynthesis_jl(temp_leaf_p::Cdouble, Rsn_leaf::Cdouble, e_air::Cdouble, - g_lb_w::Cdouble, vc_opt::Cdouble, - β_soil::Cdouble, b_h2o::Cdouble, m_h2o::Cdouble, +@fastmath function photosynthesis_jl(T_leaf_p::Cdouble, Rsn_leaf::Cdouble, ea::Cdouble, + gb_w::Cdouble, Vcmax25::Cdouble, + β_soil::Cdouble, g0_h2o::Cdouble, g1_h2o::Cdouble, cii::Cdouble, T_leaf::Cdouble, LH_leaf::Cdouble) - ca::Cdouble = CO2_air # atmospheric co2 concentration (ppm) + ca::Cdouble = CO2_air # atmospheric co2 concentration (ppm) PPFD::Cdouble = 4.55 * 0.5 * Rsn_leaf # incident photosynthetic photon flux density (PPFD) umol m-2 s-1 (2PPFD < 1) && (PPFD = 0.0) T_leaf_K = T_leaf + 273.13 - fact_latent = LAMBDA(temp_leaf_p) - bound_vapor = 1.0 / g_lb_w + λ = LAMBDA(T_leaf_p) # [J kg-1], ~2.5MJ kg-1 + bound_vapor = 1.0 / gb_w # air_pres::Cdouble = 101.325 # air pressure (kPa) - # g_lb_c = (g_lb_w/1.6)*air_pres/(temp_leaf_K*rugc); # (mol m-2 s-1) + # gb_c = (g_lb_w/1.6)*air_pres/(temp_leaf_K*rugc); # (mol m-2 s-1) atm = 1.013 # 1 atm, [1013.25 hPa] -> 1.013 bar pstat273 = 0.022624 / (273.16 * atm) - T_Kelvin = T_leaf + 273.13 - rhova_g = e_air * 2165 / T_Kelvin # absolute humidity, g m-3 - rhova_kg = rhova_g / 1000.0 # absolute humidity, kg m-3 + function cal_rho_a(Ta, ea) + TK = Ta + 273.13 + rhova_g = ea * 2165 / TK # absolute humidity, g m-3 + return rhova_g / 1000.0 # [g m-3] -> [kg m-3] + end - g_lb_c = 1.0 / (1.0 / g_lb_w * 1.6 * T_leaf_K * (pstat273)) + ρₐ = cal_rho_a(T_leaf, ea) # [kg m-3] + gb_c = 1.0 / (1.0 / gb_w * 1.6 * T_leaf_K * (pstat273)) - m_co2 = m_h2o / 1.6 # - b_co2 = b_h2o / 1.6 + g1_co2 = g1_h2o / 1.6 # + g0_co2 = g0_h2o / 1.6 - RH_leaf = SFC_VPD(T_leaf_K, LH_leaf, fact_latent, bound_vapor, rhova_kg) - tprime25 = T_leaf_K - tk_25 # temperature difference + RH_leaf = SFC_VPD(T_leaf_K, LH_leaf, λ, bound_vapor, ρₐ) + tprime25 = T_leaf_K - TK25 # temperature difference #/***** Use Arrhenius Eq. to compute KC and km_o2 *****/ - km_co2 = TEMP_FUNC(kc25, ekc, tprime25, tk_25, T_leaf_K) - km_o2 = TEMP_FUNC(ko25, eko, tprime25, tk_25, T_leaf_K) - tau = TEMP_FUNC(tau25, ektau, tprime25, tk_25, T_leaf_K) + km_co2 = TEMP_FUNC(kc25, ekc, tprime25, TK25, T_leaf_K) + km_o2 = TEMP_FUNC(ko25, eko, tprime25, TK25, T_leaf_K) + tau = TEMP_FUNC(tau25, ektau, tprime25, TK25, T_leaf_K) K = km_co2 * (1.0 + o2 / km_o2) Γ = 0.5 * o2 / tau * 1000.0 # umol mol-1 - Rd25 = vc_opt * 0.004657 #leaf dark respiration (umol m-2 s-1) + Rd25 = Vcmax25 * 0.004657 #leaf dark respiration (umol m-2 s-1) # Bin Chen: check this later. reduce respiration by 40% in light according to Amthor (2PPFD > 10) && (Rd25 *= 0.4) - Rd = TEMP_FUNC(Rd25, erd, tprime25, tk_25, T_leaf_K) + Rd = TEMP_FUNC(Rd25, erd, tprime25, TK25, T_leaf_K) - # jmopt = 29.1 + 1.64*vc_opt; Chen 1999, Eq. 7 - jmopt = 2.39 * vc_opt - 14.2 + # jmopt = 29.1 + 1.64*Vcmax25; Chen 1999, Eq. 7 + jmopt = 2.39 * Vcmax25 - 14.2 Jmax = TBOLTZ(jmopt, ejm, toptjm, T_leaf_K) # Apply temperature correction to JMAX - Vcmax = TBOLTZ(vc_opt, evc, toptvc, T_leaf_K) # Apply temperature correction to vcmax + Vcmax = TBOLTZ(Vcmax25, evc, toptvc, T_leaf_K) # Apply temperature correction to vcmax # /* - # * APHOTO = PG - resp_ld, net photosynthesis is the difference + # * APHOTO = PG - Rd, net photosynthesis is the difference # * between gross photosynthesis and dark respiration. Note # * photorespiration is already factored into PG. # # * Gs from Ball-Berry is for water vapor. It must be divided # * by the ratio of the molecular diffusivities to be valid for A # */ - α = 1.0 + (b_co2 / g_lb_c) - m_co2 * RH_leaf * β_soil - β = ca * (g_lb_c * m_co2 * RH_leaf * β_soil - 2.0 * b_co2 - g_lb_c) - γ = ca * ca * g_lb_c * b_co2 - θ = g_lb_c * m_co2 * RH_leaf * β_soil - b_co2 + α = 1.0 + (g0_co2 / gb_c) - g1_co2 * RH_leaf * β_soil + β = ca * (gb_c * g1_co2 * RH_leaf * β_soil - 2.0 * g0_co2 - gb_c) + γ = ca * ca * gb_c * g0_co2 + θ = gb_c * g1_co2 * RH_leaf * β_soil - g0_co2 # * Test for the minimum of Wc and Wj. Both have the form: # * `W = (a ci - ad)/(e ci + b)` @@ -128,20 +160,19 @@ end wc = Vcmax * (cii - Γ) / (cii + K) if (wj < wc) - # for Harley and Farquhar type model for Wj - ps_guess = wj + # A_guess = wj # Harley and Farquhar type model for Wj a_ps = Jₓ b_ps = 8.0 * Γ e_ps = 4.0 else - ps_guess = wc - a_ps = Vcmax + # A_guess = wc + a_ps = Vcmax # x1 b_ps = K e_ps = 1.0 end d_ps = Γ - # If `wj` or `wc` are less than resp_ld then A would probably be less than 0. + # If `wj` or `wc` are less than Rd then A would probably be less than 0. # This would yield a negative stomatal conductance. In this case, assume `gs` # equals the cuticular value. This assumptions yields a quadratic rather than # cubic solution for A. @@ -182,7 +213,7 @@ end # net photosynthesis rate limited by sucrose synthesis (umol m-2 s-1) j_sucrose = Vcmax / 2.0 - Rd An = min(An, j_sucrose) - + # Forests are hypostomatous. Hence, we don't divide the total resistance by 2 # since transfer is going on only one side of a leaf. @@ -198,16 +229,17 @@ end # if gs = ax + b. Use quadratic case when A <=0 # # Bin Chen: - # r_tot = 1.0/b_co2 + 1.0/g_lb_c; # total resistance to CO2 (m2 s mol-1) - # denom = g_lb_c * b_co2; + # r_tot = 1.0/b_co2 + 1.0/gb_c; # total resistance to CO2 (m2 s mol-1) + # denom = gb_c * b_co2; + # # Aquad = r_tot * e_ps; - # Bquad = (e_ps*resp_ld + a_ps)*r_tot - b_ps - e_ps*ca; - # Cquad = a_ps*(ca-d_ps) - resp_ld*(e_ps*ca+b_ps); + # Bquad = (e_ps*Rd + a_ps)*r_tot - b_ps - e_ps*ca; + # Cquad = a_ps*(ca-d_ps) - Rd*(e_ps*ca+b_ps); # original version @label quad - ps_1 = ca * g_lb_c * b_co2 - delta_1 = b_co2 + g_lb_c - denom = g_lb_c * b_co2 + ps_1 = ca * gb_c * g0_co2 + delta_1 = g0_co2 + gb_c + denom = gb_c * g0_co2 Aquad = delta_1 * e_ps Bquad = -ps_1 * e_ps - a_ps * delta_1 + e_ps * Rd * delta_1 - b_ps * denom @@ -221,9 +253,9 @@ end @label OUTDAT An = max(0.0, An) - cs = ca - An / g_lb_c + cs = ca - An / gb_c - gs_h2o_mole = (β_soil * m_h2o * RH_leaf * An / cs) + b_h2o # mol m-2 s-1 + gs_h2o_mole = (β_soil * g1_h2o * RH_leaf * An / cs) + g0_h2o # mol m-2 s-1 gs_co2_mole = gs_h2o_mole / 1.6 ci = cs - An / gs_co2_mole @@ -234,12 +266,13 @@ end + @fastmath function SFC_VPD(T_leaf_K::Float64, leleafpt::Float64, - fact_latent::Float64, bound_vapor::Float64, rhova_kg::Float64)::Float64 - es = ES(T_leaf_K) - ρv = (leleafpt / (fact_latent)) * bound_vapor + rhova_kg # kg m-3 - e = ρv * T_leaf_K / 0.2165 # mb - vpd = es - e # mb + λ::Float64, bound_vapor::Float64, ρₐ::Float64)::Float64 + es = ES(T_leaf_K) # mb + ρᵥ = (leleafpt / (λ)) * bound_vapor + ρₐ # kg m-3 + e = ρᵥ * T_leaf_K / 0.2165 # mb + vpd = es - e # mb rh = 1.0 - vpd / es # 0 to 1.0 return rh end @@ -249,10 +282,10 @@ end end -function LAMBDA(tak::Float64)::Float64 - y = 3149000.0 - 2370.0 * tak +function LAMBDA(TK::Float64)::Float64 + y = 3149000.0 - 2370.0 * TK # J kg-1 # add heat of fusion for melting ice - if tak < 273.0 + if TK < 273.0 y += 333.0 end return y