Skip to content
Snippets Groups Projects
Commit 9097a751 authored by HERIG COIMBRA Pedro Henrique's avatar HERIG COIMBRA Pedro Henrique
Browse files

xlsx + simplify imports

parent bcc09743
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
.gitignore
+ 1
− 2
View file @ 9097a751
data/tmp
data/sat
data/*/
__pycache__
scripts/__pycache__
scripts/FFP_Python/__pycache__
......
.png 0 → 100644
+ 0
− 0
View file @ 9097a751
File added
README.md
+ 13
− 3
View file @ 9097a751
......@@ -2,9 +2,19 @@
## Getting started
> open Anaconda Prompt<br>
> run 'cd PATH_TO_GIT_FOLDER'<br>
> run 'python run.py'
> cd PATH_TO_GIT_FOLDER
A. open Anaconda Prompt<br>
> conda create -n interface_sat_for_icos python=3.8<br>
> conda activate interface_sat_for_icos
B. open cmd<br>
> pip install virtualenv<br>
> virtualenv interface_sat_for_icos<br>
> source interface_sat_for_icos/bin/activate
> pip install -r requirements.txt<br>
> python run.py
## Downloading satelite data
......
data/landuse.xlsx
+ 0
− 0
View file @ 9097a751
No preview for this file type
data/stations.xlsx
+ 0
− 0
View file @ 9097a751
No preview for this file type
requirements.txt 0 → 100644
+ 39
− 0
View file @ 9097a751
affine==2.3.1
attrs==22.1.0
branca==0.6.0
certifi @ file:///C:/b/abs_ac29jvt43w/croot/certifi_1665076682579/work/certifi
charset-normalizer==2.1.1
click==8.1.3
click-plugins==1.1.1
cligj==0.7.2
colorama==0.4.6
contourpy==1.0.6
cycler==0.11.0
et-xmlfile==1.1.0
Fiona==1.8.22
folium==0.13.0
fonttools==4.38.0
geopandas==0.12.1
idna==3.4
Jinja2==3.1.2
kiwisolver==1.4.4
MarkupSafe==2.1.1
matplotlib==3.6.2
munch==2.5.0
numpy==1.23.4
openpyxl==3.0.10
packaging==21.3
pandas==1.5.1
Pillow==9.3.0
pyparsing==3.0.9
pyproj==3.4.0
python-dateutil==2.8.2
pytz==2022.6
rasterio==1.3.4
requests==2.28.1
scipy==1.9.3
Shapely==1.8.5.post1
six==1.16.0
snuggs==1.4.7
urllib3==1.26.12
wincertstore==0.2
run.py
+ 34
− 22
View file @ 9097a751
......@@ -2,8 +2,10 @@ import webbrowser
import scripts.common as tt
import scripts.satelite_products as satprod
import scripts.satelite as satelite
import pandas as pd
import numpy as np
from pandas import to_datetime, DataFrame, read_excel, unique
from numpy import nan, abs as npabs
from numpy.random import random as nprand, normal as npnorm
from numpy.ma import array as maarray
import tkinter as tk
import tkinter.ttk as ttk
......@@ -131,7 +133,7 @@ def plot_map(widget=None, menu_options=[]):
if widget is None:
widget = currentPreview
datetime2plot = datetime2plot if datetime2plot else pd.to_datetime('20210301')
datetime2plot = datetime2plot if datetime2plot else to_datetime('20210301')
map2plot = 'LAND USE' if not rmap2plot.get() else rmap2plot.get().upper()
site_ec = tt.ECsite(listSitelb.get(tk.ANCHOR))
......@@ -171,18 +173,19 @@ def plot_map(widget=None, menu_options=[]):
bgmap = satprod.satprod().get(map2plot, center=site_ec.loc, buffer=max(bufferMap.get(), 200),
timelist=datetime2plot, tile=site_ec.tile, download=False, filepath=site_ec.lu_path,
verbosity=0, **vprm_needs,
satfolder=wrkpath.get("1.0", "end-1c") + '/sat/SENTINEL2/')
satfolder=satpath.get("1.0", "end-1c") + '/SENTINEL2/')
if FPactivate.get()==1:
from scripts.footprint import footprint
FP = footprint(loc=site_ec.loc,
df=pd.DataFrame({'zm': [FPinput_Zm.get()],
'z0': [FPinput_Zm.get()/1000],
'h': [1000],
'ol': [FPinput_OL.get()],
'sigmav': [FPinput_SV.get()],
'ustar': [FPinput_US.get()],
'wd': [FPinput_WD.get()]}))
df=DataFrame({'zm': [FPinput_Zm.get()]*FPinput_N.get(),
'z0': [FPinput_Zm.get()/1000]*FPinput_N.get(),
'h': [1000]*FPinput_N.get(),
'ol': [FPinput_OL.get()]*FPinput_N.get(),
'sigmav': [FPinput_SV.get()]*FPinput_N.get(),
'ustar': [FPinput_US.get()] if FPinput_N.get() == 1 else list(npabs(npnorm(FPinput_US.get(), 0.1, FPinput_N.get()))),
'wd': [FPinput_WD.get()] if FPinput_N.get() == 1 else list(nprand(FPinput_N.get())*360)
}))
FP.get(r=[10,20,30,40,50,60,70,80,90], dx=bgmap.dx, dy=bgmap.dy,
nx=bgmap.data.shape[-1]-1, ny=bgmap.data.shape[-2]-1, verbosity=0)
FP = FP.FP
......@@ -196,7 +199,7 @@ def plot_map(widget=None, menu_options=[]):
lu_data = satprod.satprod().get('LU', center=site_ec.loc, buffer=max(bufferMap.get(), 200),
filepath=site_ec.lu_path, **kwargs).data
lu_meta = pd.read_excel(site_ec.lu_metapath).set_index(
lu_meta = read_excel(site_ec.lu_metapath).set_index(
'code')[['name', 'IGBP']].to_dict()
lu_dic = {}
for k, v in lu_meta['IGBP'].items(): # <- missing ()
......@@ -204,9 +207,9 @@ def plot_map(widget=None, menu_options=[]):
lu_dic[v].append(k)
lu_data = satprod.force_array_dimension(bgmap.data[0].shape, lu_data)
lu_data = sum([np.ma.array(lu_data, mask=lu_data == el)
lu_data = sum([maarray(lu_data, mask=lu_data == el)
for el in lu_dic[lu_]])
lu_data[~lu_data.mask] = np.nan
lu_data[~lu_data.mask] = nan
lu_data[lu_data.mask] = 1
lu_filter = tt.sum_nan_arrays(lu_filter, lu_data) if 'lu_filter' in locals() else lu_data
......@@ -260,7 +263,7 @@ def luChoice(lu):
def dateChoice(e=None):
global datetime2plot
try:
datetime2plot = pd.to_datetime(listDatelb.get(tk.ANCHOR))
datetime2plot = to_datetime(listDatelb.get(tk.ANCHOR))
#btdate2plot.configure(text='-> ' + str(datetime2plot))
except ValueError:
datetime2plot = None
......@@ -273,12 +276,12 @@ def updateDateList(e=None):
listDatelb.delete(0, tk.END)
try:
[listDatelb.insert(i, el) for i, el in enumerate(sorted(set(satelite.get_nearest_satelite(
wrkpath.get("1.0", "end-1c") + '/sat/SENTINEL2/', tile=tt.ECsite(listSitelb.get(tk.ANCHOR)).tile)[1])))]
satpath.get("1.0", "end-1c") + '/SENTINEL2/', tile=tt.ECsite(listSitelb.get(tk.ANCHOR)).tile)[1])))]
except:
return
def countryChoice(e=None):
site_selection = pd.read_excel(wrkpath.get("1.0", "end-1c")+"/stations.xlsx")
site_selection = read_excel(wrkpath.get("1.0", "end-1c")+"/stations.xlsx")
if countryChosen.cget("text"):
site_selection = site_selection[site_selection.CO == countryChosen.cget("text")]
site_selection = tt.prioritize_list(site_selection.Name.to_list(), "(Grig)|(Hegyhátsál)|(Sac)|(Ris)")
......@@ -295,7 +298,7 @@ def updateSiteList(e=None):
"""
Globals
"""
datetime2plot = pd.to_datetime('20210301')
datetime2plot = to_datetime('20210301')
root = tk.Tk()
root.grid_columnconfigure(1, weight=1)
......@@ -361,6 +364,12 @@ wrkpath = tk.Text(display_wrk, height=1, width=50)
wrkpath.insert(tk.END, 'data/') # ../eddy-tower_largefiles/data/input/raw
wrkpath.grid(row=0, column=1, sticky=tk.W)
tk.Label(display_wrk, text='Sat. folder:').grid(row=1, column=0)
satpath = tk.Text(display_wrk, height=1, width=50)
satpath.insert(tk.END, 'data/sat/') # ../eddy-tower_largefiles/data/input/raw
satpath.grid(row=1, column=1, sticky=tk.W)
"""
Main Menu ___________
"""
......@@ -403,7 +412,7 @@ listCountrylb.grid(row=0, column=0, sticky=tk.W+tk.E)
listCountrylb.configure(yscrollcommand=listCountrysb.set)
listCountrysb.config(command=listCountrylb.yview)
[listCountrylb.insert(i, el) for i, el in enumerate([''] + tt.prioritize_list(list(pd.unique(pd.read_excel(
[listCountrylb.insert(i, el) for i, el in enumerate([''] + tt.prioritize_list(list(unique(read_excel(
wrkpath.get("1.0", "end-1c")+"/stations.xlsx").CO)), "(FR)|(DK)|(HU)"))]
listCountrylb.bind('<<ListboxSelect>>', lambda e: countryChosen.configure(
......@@ -444,7 +453,7 @@ listSitelb.bind('<<ListboxSelect>>', updateDateList, '+')
listSitelb.select_anchor(0)
listDatelb.bind('<<ListboxSelect>>', dateChoice)
listDatelb.bind('<<ListboxSelect>>', lambda e: dateChosen.configure(text=pd.to_datetime(e.widget.get(
listDatelb.bind('<<ListboxSelect>>', lambda e: dateChosen.configure(text=to_datetime(e.widget.get(
e.widget.curselection())).date()) if e.widget.curselection() else None, '+')
......@@ -452,7 +461,7 @@ def winDownloadSatDate():
new_window()
nwind.title('VPRM Equations')
tile = pd.read_excel(
tile = read_excel(
wrkpath.get("1.0", "end-1c")+"/stations.xlsx")
tk.Label(nwind, text=listSitelb.get(tk.ANCHOR), justify='left').pack()
tile = tile[tile.Name == listSitelb.get(tk.ANCHOR)].SateliteTile.values[0]
......@@ -607,7 +616,10 @@ innerframeFP.grid_columnconfigure(1, weight=1)
"""
FPactivate = tk.IntVar()
tk.Checkbutton(innerframeFP, text="Calculate footprint.", variable=FPactivate,
onvalue = 1, offvalue = 0).grid(row=0, column=1, sticky=tk.W)
onvalue=1, offvalue=0).grid(row=0, column=1, sticky=tk.W)
FPinput_N = tk.Scale(innerframeFP, from_=1, to=48,
orient=tk.HORIZONTAL)
FPinput_N.grid(row=0, column=1, sticky=tk.E)
tk.Label(innerframeFP, text='Zm (m): ').grid(row=1, column=0)
FPinput_Zm = tk.Scale(innerframeFP, from_=2, to=120,
......
scripts/common.py
+ 84
− 192
View file @ 9097a751
This diff is collapsed.
scripts/satelite.py
+ 19
− 18
View file @ 9097a751
import scripts.common as tcom
import pandas as pd
from rasterio.io import MemoryFile
from rasterio import mask as msk
from rasterio import plot as rastplot
import rasterio
from rasterio import mask as msk, plot as rastplot, open as rastopen, band
from rasterio.warp import reproject
from shapely import geometry, ops
from pyproj import Transformer
import os
import re
import zipfile
import numpy as np
from numpy import isnan, nan, sin, cos, radians, zeros
from matplotlib import cm, pyplot
from matplotlib import cm, pyplot, lines as mplLine
from matplotlib.colors import rgb2hex
import matplotlib.lines as mplLine
from fiona.drvsupport import supported_drivers
supported_drivers['KML'] = 'rw'
#def download_satdata(stdate='2021-03-02', endate='2021-03-31', tile='T31UDQ', sat='SENTINEL2', cloudcover='50'):
......@@ -34,11 +35,11 @@ def download_satdata(stdate, endate, tile, sat='SENTINEL2', cloudcover='50'):
def distfromtower(center, countour, wd, verbosity=0):
if np.isnan(wd):
return np.nan
if isnan(wd):
return nan
ln = geometry.LineString(
[(center), (center[0] + 500 * np.sin(np.radians(wd)),
center[1] + 500 * np.cos(np.radians(wd)))])
[(center), (center[0] + 500 * sin(radians(wd)),
center[1] + 500 * cos(radians(wd)))])
intersec = countour.intersection(ln)
if verbosity:
pyplot.plot(*ln.xy)
......@@ -50,8 +51,8 @@ def distfromtower(center, countour, wd, verbosity=0):
def crop_tif_by_polygon(cutshape, basegrid, extent):
# time consuming !
# crop meshgrid by checking each point if it is inside a polygon
outputgrid = np.zeros(basegrid.shape) #.astype(float)
outputgrid[:] = np.nan
outputgrid = zeros(basegrid.shape) #.astype(float)
outputgrid[:] = nan
for i in range(len(basegrid)):
for j in range(len(basegrid[i])):
......@@ -94,7 +95,7 @@ def cut_proj(infile, center, buffer, extentBox=None, outfolder="data/tmp/", outf
[center[0]+buffer, center[1]+buffer],
[center[0]-buffer, center[1]+buffer]])
with rasterio.open(infile) as src:
with rastopen(infile) as src:
# reproject slighlty bigger box
transformer = Transformer.from_crs(
"EPSG:3035", src.crs, always_xy=True).transform
......@@ -109,8 +110,8 @@ def cut_proj(infile, center, buffer, extentBox=None, outfolder="data/tmp/", outf
src_cut.write(out_cut_image)
# reproject to 3035
src_reproj, src_reproj_trans = rasterio.warp.reproject(
source=rasterio.band(src_cut, 1), dst_crs='EPSG:3035')
src_reproj, src_reproj_trans = reproject(
source=band(src_cut, 1), dst_crs='EPSG:3035')
#src_reproj = src_reproj[0]
# translate back to tif format
......@@ -129,7 +130,7 @@ def cut_proj(infile, center, buffer, extentBox=None, outfolder="data/tmp/", outf
# , "crs": pycrs.parser.from_epsg_code(epsg_code).to_proj4()})
"""
with rasterio.open(infile) as src:
with rastopen(infile) as src:
transformer = Transformer.from_crs("EPSG:3035", src.crs, always_xy=True).transform
gbox = ops.transform(transformer, gbox)
# print(gbox.bounds)
......@@ -139,12 +140,12 @@ def cut_proj(infile, center, buffer, extentBox=None, outfolder="data/tmp/", outf
#print(out_meta)
"""
with rasterio.open(outfile, "w", **out_meta) as dest:
with rastopen(outfile, "w", **out_meta) as dest:
dest.write(out_image)
cutproj_meta.write()
clipped = rasterio.open(outfile)
clipped = rastopen(outfile)
if plot:
rastplot.show((clipped, 1), cmap='terrain')
......
scripts/satelite_products.py
+ 21
− 41
View file @ 9097a751
import numpy as np
import scripts.common as tt
from scripts.satelite import *
from scripts.surfmodel import vprm, vprm_
import scipy as sp
from pandas import to_datetime
from scripts.common import sum_nan_arrays, metadata
from scripts.satelite import cut_proj, crop_tif_by_polygon, plot_sat, get_path_to_satelite, get_sat_data
from rasterio import open as rastopen
from rasterio.plot import plotting_extent
from scripts.surfmodel import vprm
from scipy.ndimage import zoom
from pandas import to_datetime, Timedelta
from os import path
"""Meta for Satelite Products"""
satproducts_dic = {
......@@ -55,7 +58,7 @@ class satprod:
"""
landuse = cut_proj(filepath, **kwargs)
self.data = landuse.read(1)
self.extent = rastplot.plotting_extent(landuse)
self.extent = plotting_extent(landuse)
self.profile = landuse.profile
elif product.upper() in ["NDVI", "EVI", "LSWI", "RGB"]:
......@@ -253,7 +256,7 @@ def force_array_dimension(shape, out):
out_[bas_cut[0]:bas_cut[1],
bas_cut[2]:bas_cut[3]] = \
tt.sum_nan_arrays(out_[bas_cut[0]:bas_cut[1],
sum_nan_arrays(out_[bas_cut[0]:bas_cut[1],
bas_cut[2]:bas_cut[3]],
out[out_cut[0]:out_cut[1],
out_cut[2]:out_cut[3]])
......@@ -263,7 +266,7 @@ def force_array_dimension(shape, out):
def __sat_product_i__(sat_foldername, sat_date, sat_date_dif, center, buffer=600, returnTIF=True, method="NDVI", **kwargs):
satprod_meta = tcom.metadata(**satproducts_dic[method])
satprod_meta = metadata(**satproducts_dic[method])
sat_data, satdata_meta = get_sat_data(sat_date, center=center, sat_foldername=sat_foldername,
sat_date_dif=sat_date_dif,
......@@ -272,15 +275,15 @@ def __sat_product_i__(sat_foldername, sat_date, sat_date_dif, center, buffer=600
band = satdata_meta.bands
satprod_meta.update(**vars(satdata_meta))
sat_date = os.path.split(sat_data[1].name)[0].split('/')[-1]
sat_date = path.split(sat_data[1].name)[0].split('/')[-1]
# Preparing metadata
prod_path = "data/tmp/SENTINEL2/" + sat_date + satprod_meta.filename
prod_metapath = os.path.splitext(prod_path)[0] + ".metadata"
prod_metapath = path.splitext(prod_path)[0] + ".metadata"
satprod_meta.filepath = prod_metapath
same_metadata = os.path.exists(prod_path) and \
same_metadata = path.exists(prod_path) and \
False #satprod_meta.check(attrs=['formula', 'datedif', 'ref', 'tile', 'buffer', 'center'])
if not same_metadata:
......@@ -307,11 +310,11 @@ def __sat_product_i__(sat_foldername, sat_date, sat_date_dif, center, buffer=600
in sentinel 2 for SWIR we have B11 and B12 both at 20 m resolution,
so we need to pass band 11 from 20 m to 10 m (1 pixel -> 4 equal pixels)
"""
shp_11 = np.array(rastplot.plotting_extent(sat_data[band['B11']]))
shp_8 = np.array(rastplot.plotting_extent(sat_data[band['B8']]))
shp_11 = np.array(plotting_extent(sat_data[band['B11']]))
shp_8 = np.array(plotting_extent(sat_data[band['B8']]))
corr_20m_to_10m = np.round((shp_8-shp_11)/10, 0).astype(int)
corr_20m_to_10m = [c if c!=0 else None for c in corr_20m_to_10m]
sat_data_B11 = sp.ndimage.zoom(sat_data[band['B11']].read(1),[2,2],order=0,mode='nearest')
sat_data_B11 = zoom(sat_data[band['B11']].read(1),[2,2],order=0,mode='nearest')
sat_data_B11 = np.array([el[corr_20m_to_10m[0]:corr_20m_to_10m[1]] for el in sat_data_B11[corr_20m_to_10m[2]:corr_20m_to_10m[3]]])
if sat_data[band['B8']].read(1).shape != sat_data_B11.shape:
......@@ -353,48 +356,25 @@ def __sat_product_i__(sat_foldername, sat_date, sat_date_dif, center, buffer=600
prod_meta.update(**{'dtype': 'float64', "nodata": np.nan})
prod_meta.update(
**{'datedif': satprod_meta.datedif/pd.Timedelta(1, 'd')})
**{'datedif': satprod_meta.datedif/Timedelta(1, 'd')})
satprod_meta.update(**prod_meta)
with rasterio.open(prod_path, "w", **prod_meta) as dest:
with rastopen(prod_path, "w", **prod_meta) as dest:
if method=="RGB": dest.write(np.array(prod_data))
else: dest.write(np.array(list([prod_data])))
satprod_meta.write()
with rasterio.open(prod_path, "r") as prod_file:
with rastopen(prod_path, "r") as prod_file:
if returnTIF:
return(prod_file)
else:
if method=="RGB": prod_data = prod_file.read()
else: prod_data = prod_file.read(1)
product = {"data": prod_data,
"extent": rastplot.plotting_extent(prod_file),
"extent": plotting_extent(prod_file),
"profile": prod_file.profile,
"meta": satprod_meta}
return product
""" ALTIMETRY
# Get cut projection
def get_altimetry(center, buffer, FP=None, plot=False):
alt_tif = cut_proj("../eddy-tower_largefiles/data/input/raw/sat/FR001L1_PARIS_UA2012_DHM_v010/Data/FR001L1_PARIS_UA2012_DHM_v010.tif",
center=center, buffer=buffer, force=False, plot=False)
# Mask nan (=-32767)
alt_arr = np.ma.masked_array(data = alt_tif.read(1),
mask = [a==-32767 for a in alt_tif.read(1)],
fill_value = -32767)
# reverse to reposition North (could be more robust)
extent = rastplot.plotting_extent(alt_tif)
extent = extent[:2] + extent[3:] + extent[2:3]
# Plot
if plot:
plot_sat(alt_arr, footprint=FP, mast=site_loc, extent=extent)
return(alt_tif)
"""
\ No newline at end of file
scripts/surfmodel.py
+ 12
− 117
View file @ 9097a751
import scripts.common as tt
import numpy as np
default_params = {
'Tmin': {'EF': 0, 'DF': 0, 'MF': 0, 'CRO': 0, 'GRA': 0},
'Topt': {'EF': 20, 'DF': 20, 'MF': 20, 'CRO': 22, 'GRA': 20},
'Tmax': {'EF': 40, 'DF': 40, 'MF': 40, 'CRO': 40, 'GRA': 40},
'Tlow': {'EF': 5, 'DF': 5, 'MF': 5, 'CRO': 2, 'GRA': 5},
'PAR0': {'EF': 275, 'DF': 254, 'MF': 446, 'CRO': 1132, 'GRA': 528},
'λ': {'EF': 0.226, 'DF': 0.215, 'MF': 0.163, 'CRO': 0.086, 'GRA': 0.119},
'α': {'EF': 0.288, 'DF': 0.181, 'MF': 0.244, 'CRO': 0.092, 'GRA': 0.125},
'β': {'EF': -1.100, 'DF': 0.840, 'MF': -0.490, 'CRO': 0.290, 'GRA': 0.017}
}
def pvprmd():
"""Get default parameters"""
result = tt.parameter_model(['Tmin', 'Topt', 'Tmax', 'Tlow', 'PAR0', 'λ', 'α', 'β'],
['EF', 'DF', 'MF', 'CRO', 'GRA'],
**{
'Tmin': {'EF': 0, 'DF': 0, 'MF': 0, 'CRO': 0, 'GRA': 0},
'Topt': {'EF': 20, 'DF': 20, 'MF': 20, 'CRO': 22, 'GRA': 20},
'Tmax': {'EF': 40, 'DF': 40, 'MF': 40, 'CRO': 40, 'GRA': 40},
'Tlow': {'EF': 5, 'DF': 5, 'MF': 5, 'CRO': 2, 'GRA': 5},
'PAR0': {'EF': 275, 'DF': 254, 'MF': 446, 'CRO': 1132, 'GRA': 528},
'λ': {'EF': 0.226, 'DF': 0.215, 'MF': 0.163, 'CRO': 0.086, 'GRA': 0.119},
'α': {'EF': 0.288, 'DF': 0.181, 'MF': 0.244, 'CRO': 0.092, 'GRA': 0.125},
'β': {'EF': -1.100, 'DF': 0.840, 'MF': -0.490, 'CRO': 0.290, 'GRA': 0.017}
})
return result
from numpy import expand_dims, array, nanmax, exp
def vprm(EVI, LSWI, PAR, T, PAR0, Tlow, Tmax, Tmin, Topt, α, β, λ, lu_map=None, lu_dic=None, verbosity=0, **kwargs):
"""
......@@ -37,25 +6,25 @@ def vprm(EVI, LSWI, PAR, T, PAR0, Tlow, Tmax, Tmin, Topt, α, β, λ, lu_map=Non
Xiao et al., 2004;
Yuan et al., 2015
"""
T = np.expand_dims(np.array(T), axis=(1, 2))
PAR = np.expand_dims(np.array(PAR), axis=(1, 2))
Tmin = np.array([Tmin]*len(T))
Tmax = np.array([Tmax]*len(T))
Topt = np.array([Topt]*len(T))
#λ = np.array([λ]*len(T))
#PAR0 = np.array([PAR0]*len(T))
#α = np.array([α]*len(T))
#β = np.array([β]*len(T))
T = expand_dims(array(T), axis=(1, 2))
PAR = expand_dims(array(PAR), axis=(1, 2))
Tmin = array([Tmin]*len(T))
Tmax = array([Tmax]*len(T))
Topt = array([Topt]*len(T))
#λ = array([λ]*len(T))
#PAR0 = array([PAR0]*len(T))
#α = array([α]*len(T))
#β = array([β]*len(T))
Tscale = ((T - Tmin) * (T - Tmax)) / \
(((T - Tmin) * (T - Tmax)) - (T - Topt)**2)
# Pscalar = (1 + LSWI)/2 [During bud burst to leaf full expansion] and = 1 [After leaf full expansion]
Pscale = (1 + LSWI) / 2
Wscale = (1 + LSWI) / (1 + np.nanmax(LSWI))
Wscale = (1 + LSWI) / (1 + nanmax(LSWI))
fgpp_ls = - (λ * Tscale * Pscale * Wscale) * \
(1 / (1 + PAR/PAR0)) * EVI * PAR
reco_exp_t = np.exp(4.6 * (1/40 - 1/T))
reco_exp_t = exp(4.6 * (1/40 - 1/T))
reco_exp_t[T[:] < 0] = 0
reco_ls_a = α * T #* reco_exp_t
reco_ls_b = β
......@@ -70,77 +39,3 @@ def vprm_(*args, **kwargs):
def save_meshgrid_tif(dat, profile):
# save it so that I can plot histogram
return
"""
def fco2_px(mask, fco2, inhomogeneity, PAR=None, hour=None, param=[1,1,0,0]):
#assert(PAR or hour)
k, r, r0, k0 = param
fco2_ls = []
photosyn = [p/np.nanmax(PAR) if p>100 else 0 for p in PAR] # the higher the PAR the more ndvi affects co2 spatial variability
for i in range(len(fco2)):
fco2_px = mask * k * inhomogeneity[i] ** (photosyn[i] + r0) * fco2[i] + k0
fco2_ls = fco2_ls + [fco2_px]
return fco2_ls
def fco2_lo(param):
T_param, PAR0, λ, Reco_param = vprm_default_values()
param = [T_param] + [PAR0] + list(param) + [Reco_param]
#_estmat = fco2_px(lumask['data'], np.array(obs)[TFlist], np.array(ndvi_mat)[TFlist], np.array(PPFD)[TFlist], np.array(hour)[TFlist], param)
_estmat = vprm(np.array(evi_mat)[TFlist], np.array(lswi_mat)[TFlist], np.array(PPFD)[TFlist],
(dft[TFlist].air_temperature-272.15).tolist(), T_param=param[0:4], PAR0=param[4], λ=param[5], Reco_param=param[6:8])
_est = np.array([np.nansum(_estmat[i]*fp_mat[i]) for i in range(len(_estmat))])
#print(np.array(_est[:time_length]).shape)
return _est #_est[:time_length]
#def vprm(EVI, LSWI, PAR, T, T_param=[0,20,40,5], PAR0=570, λ=0.127, Reco_param=[0.271, 0.25])
def fco2_lo_r(obs_):
return
def fco2_lo_r(dat, evi_ls, lswi_ls, par_ls, ta_ls, fp_ls):
def compute_distance(param, data):
# parameters
#T_param, PAR0, λ, Reco_param = vprm_default_values()
#param = [T_param] + [PAR0] + list(param) + [Reco_param]
param = vprm_default_values(λ=param)
EVI, LSWI, PAR, T, fco2, footprint = data
# run vprm
estmat = vprm(EVI, LSWI, PAR, T, T_param=param[0], PAR0=param[1], λ=param[2], Reco_param=param[3])
# turn meshgrid into timeseries (w/ footprint)
est = np.array([np.nansum(estmat[i]*footprint[i]) for i in range(len(estmat))])
# get diff we try to minimize
distance = np.sum((fco2-est)**2)
return(distance)
bnds = ((0, 1),)
opt = sp.optimize.minimize(compute_distance, [vprm_default_values()[2]],
args=([evi_ls, lswi_ls, par_ls, ta_ls, dat, fp_ls],),
bounds=bnds)
return opt.x
def fco2_lo_rBKP(dat, ):
#print("s", end='')
#bnds = ((0, 1.5), (0, 1.5), (-1.5, 1.5), (-1.5, 1.5))
bnds = ((0, 1))
opt = sp.optimize.minimize(compute_distance,
[vprm_default_values()[2]],# param_opt,
#args=([lumask['data'], np.array(obs_), np.array(ndvi_mat)[TFlist], np.array(PPFD)[TFlist], np.array(hour)[TFlist], np.array(fp_mat)[TFlist]],),
args=([np.array(evi_mat)[TFlist], np.array(lswi_mat)[TFlist], np.array(PPFD)[TFlist],
(dft[TFlist].air_temperature-272.15).tolist(), np.array(obs_), np.array(fp_mat)[TFlist]],),
#bounds=bnds
)
#print(opt.x)
#print(param)
#_est = fco2_lo(param)
#_est_r = np.array(list(map(list, zip(*_est))))
#print("o")
return opt.x
"""
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment