Skip to content

use fiona-based antimeridian-safe crop and reprojection for land cover map #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Timothy-W-Hilton wants to merge 4 commits into
base: main
Choose a base branch
from
Open

use fiona-based antimeridian-safe crop and reprojection for land cover map #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
514a494
blacken python source files
Timothy-W-Hilton Sep 23, 2024
727602a
use loguru package for logging
Timothy-W-Hilton Sep 23, 2024
5aa2fd7
BUGFIX: use cartesian coordinates for disjoint bounds test
Timothy-W-Hilton Sep 24, 2024
692608c
use fiona-based antimeridian-safe crop for LCM
Timothy-W-Hilton Oct 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 47 additions & 30 deletions fit_vprm_parameters/fit_params_draft.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,42 +1,44 @@
import sys
import os
import pathlib
sys.path.append(os.path.join(pathlib.Path(__file__).parent.resolve(), '..', 'lib'))

sys.path.append(os.path.join(pathlib.Path(__file__).parent.resolve(), "..", "lib"))
from sat_manager import VIIRS, modis, copernicus_land_cover_map
from VPRM import vprm
from VPRM import vprm
import yaml
import glob
import time
import numpy as np
import pandas as pd
import pickle
import argparse
from loguru import logger

p = argparse.ArgumentParser(
description = "Commend Line Arguments",
formatter_class = argparse.RawTextHelpFormatter)
description="Commend Line Arguments", formatter_class=argparse.RawTextHelpFormatter
)
p.add_argument("--h", type=int)
p.add_argument("--v", type=int)
p.add_argument("--config", type=str)
p.add_argument("--n_cpus", type=int, default=1)
args = p.parse_args()
print(args)
logger.info(str(args))


with open(args.config, "r") as stream:
try:
cfg = yaml.safe_load(stream)
cfg = yaml.safe_load(stream)
except yaml.YAMLError as exc:
print(exc)
logger.info(exc)

h = args.h
v = args.v

if not os.path.exists(cfg['out_path']):
os.makedirs(cfg['out_path'])
if not os.path.exists(cfg["out_path"]):
os.makedirs(cfg["out_path"])

outfile = os.path.join(cfg['out_path'], 'h{}v{}_{}.pickle'.format(h, v, cfg['year']))
print(outfile)
outfile = os.path.join(cfg["out_path"], "h{}v{}_{}.pickle".format(h, v, cfg["year"]))
logger.info(outfile)


### ToDo: Provide a list of lats and lons for the data extraction, i.e. the flux tower positions
Expand All @@ -51,25 +53,37 @@

vprm_inst = vprm(n_cpus=args.n_cpus)

for c, i in enumerate(glob.glob(os.path.join(cfg['sat_image_path'], '*h{:02d}v{:02d}*.h*'.format(h, v)))):
print(i)
if cfg['satellite'] == 'modis':
for c, i in enumerate(
glob.glob(os.path.join(cfg["sat_image_path"], "*h{:02d}v{:02d}*.h*".format(h, v)))
):
logger.info(i)
if cfg["satellite"] == "modis":
handler = modis(sat_image_path=i)
handler.load()
vprm_inst.add_sat_img(handler, b_nir='B02', b_red='B01',
b_blue='B03', b_swir='B06',
which_evi='evi',
drop_bands=True)
vprm_inst.add_sat_img(
handler,
b_nir="B02",
b_red="B01",
b_blue="B03",
b_swir="B06",
which_evi="evi",
drop_bands=True,
)
else:
handler = VIIRS(sat_image_path=i)
handler.load()
vprm_inst.add_sat_img(handler, b_nir='SurfReflect_I2', b_red='SurfReflect_I1',
b_blue='no_blue_sensor', b_swir='SurfReflect_I3',
which_evi='evi2',
drop_bands=True)
vprm_inst.add_sat_img(
handler,
b_nir="SurfReflect_I2",
b_red="SurfReflect_I1",
b_blue="no_blue_sensor",
b_swir="SurfReflect_I3",
which_evi="evi2",
drop_bands=True,
)

# Uncomment if you want to do a spatial smearing with size=(3,3). Otherwise give size manually.
#vprm_inst.smearing(lonlats=lonlats)
# vprm_inst.smearing(lonlats=lonlats)

# Sort the satellite images by time and merge internally for easier computations
vprm_inst.sort_and_merge_by_timestamp()
Expand All @@ -84,15 +98,18 @@

# Loop over the year and get the vprm_variables

time_range = pd.date_range(start="{}-01-01".format(cfg['year']),
end="{}-01-01".format(cfg['year'] + 1),
freq='H')
time_range = pd.date_range(
start="{}-01-01".format(cfg["year"]),
end="{}-01-01".format(cfg["year"] + 1),
freq="H",
)

for c,t in enumerate(time_range):
for c, t in enumerate(time_range):
t0 = time.time()
# Returns dict with VPRM input variables for the different lats and lons
# Additional ERA5 variables can be added trough the add_era_variables key
vrbls = vprm_inst.get_vprm_variables(t, lat=lats, lon=lons,
add_era_variables=['svwl1'])

vrbls = vprm_inst.get_vprm_variables(
t, lat=lats, lon=lons, add_era_variables=["svwl1"]
)

# ToDo: write to your own format and store for further analysis
150 changes: 96 additions & 54 deletions fit_vprm_parameters/get_fit_files_for_site.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
from pyVPRM.lib.flux_tower_class import fluxnet, icos
from pyVPRM.sat_managers.viirs import VIIRS
from pyVPRM.sat_managers.modis import modis
from pyVPRM.VPRM import vprm
from pyVPRM.VPRM import vprm
import xarray as xr
import pickle
import yaml
Expand All @@ -15,6 +15,8 @@
import datetime
from astropy.convolution import Gaussian2DKernel
from functions import lat_lon_to_modis
from loguru import logger


def all_files_exist(item):
for key in item.assets.keys():
Expand All @@ -23,9 +25,10 @@ def all_files_exist(item):
return False
return True


p = argparse.ArgumentParser(
description = "Commend Line Arguments",
formatter_class = argparse.RawTextHelpFormatter)
description="Commend Line Arguments", formatter_class=argparse.RawTextHelpFormatter
)
p.add_argument("--site", type=str)
p.add_argument("--veg_type", type=str)
p.add_argument("--this_year", type=str)
Expand All @@ -34,100 +37,137 @@ def all_files_exist(item):

with open(args.cfg_path, "r") as stream:
try:
cfg = yaml.safe_load(stream)
cfg = yaml.safe_load(stream)
except yaml.YAMLError as exc:
print(exc)
logger.info(exc)

s = args.site
veg_type = args.veg_type
this_year = int(args.this_year)
this_year = int(args.this_year)

veg_type_id = {'GRA': 7, 'MF': 3 , 'CRO':6,
'EF': 1, 'DF': 2, 'SH':4, 'WET': 9}
veg_type_id = {"GRA": 7, "MF": 3, "CRO": 6, "EF": 1, "DF": 2, "SH": 4, "WET": 9}

# GRA = Grassland | MF = Mixed Forest | CRO = Cropland | EF = Evergreen Forest | DF = Deciduous Forest | SH = Shrubland

all_data = []
vprm_insts = []
tmin = parser.parse('{}0101'.format(this_year))
tmax = parser.parse('{}1231'.format(this_year))
tmin = parser.parse("{}0101".format(this_year))
tmax = parser.parse("{}1231".format(this_year))
tower_data_list = []
print(s)
logger.info(s)
# if len(tower_data_list) >1:
# continue
if this_year == 2012:
data_files = glob.glob(os.path.join(cfg['fluxnet_path'], '*{}*/*_FULLSET_H*'.format(s)))
if len(data_files) == 0: exit()
flux_tower_inst = fluxnet(data_files[0],'SW_IN_F', 'TA_F',
t_start= tmin, t_stop=tmax)
data_files = glob.glob(
os.path.join(cfg["fluxnet_path"], "*{}*/*_FULLSET_H*".format(s))
)
if len(data_files) == 0:
exit()
flux_tower_inst = fluxnet(
data_files[0], "SW_IN_F", "TA_F", t_start=tmin, t_stop=tmax
)
elif this_year in [2021, 2022]:
data_files = glob.glob(os.path.join(cfg['icos_path'], '*/*_{}_FLUXNET_HH_L2.csv'.format(s)))
if len(data_files) == 0: exit()
flux_tower_inst = icos(data_files[0],'SW_IN_F', 'TA_F', t_start= tmin, t_stop=tmax)
data_files = glob.glob(
os.path.join(cfg["icos_path"], "*/*_{}_FLUXNET_HH_L2.csv".format(s))
)
if len(data_files) == 0:
exit()
flux_tower_inst = icos(data_files[0], "SW_IN_F", "TA_F", t_start=tmin, t_stop=tmax)
lon = flux_tower_inst.get_lonlat()[0]
lat= flux_tower_inst.get_lonlat()[1]
lat = flux_tower_inst.get_lonlat()[1]
okay = flux_tower_inst.add_tower_data()
if not okay: exit()
print(flux_tower_inst.get_site_name())
if not okay:
exit()
logger.info(flux_tower_inst.get_site_name())
flux_tower_inst.set_land_type(veg_type_id[veg_type])

lat = flux_tower_inst.get_lonlat()[1]
lon = flux_tower_inst.get_lonlat()[0]
h, v = lat_lon_to_modis(lat, lon)
inp_files1 = sorted(glob.glob(os.path.join(cfg['sat_image_path'], str(this_year-1),
'*h{:02d}v{:02d}*.h*'.format(h, v))))[-8:]
inp_files2 = sorted(glob.glob(os.path.join(cfg['sat_image_path'], str(this_year),
'*h{:02d}v{:02d}*.h*'.format(h, v))))
inp_files3 = sorted(glob.glob(os.path.join(cfg['sat_image_path'], str(this_year+1),
'*h{:02d}v{:02d}*.h*'.format(h, v))))[:8]
inp_files1 = sorted(
glob.glob(
os.path.join(
cfg["sat_image_path"],
str(this_year - 1),
"*h{:02d}v{:02d}*.h*".format(h, v),
)
)
)[-8:]
inp_files2 = sorted(
glob.glob(
os.path.join(
cfg["sat_image_path"], str(this_year), "*h{:02d}v{:02d}*.h*".format(h, v)
)
)
)
inp_files3 = sorted(
glob.glob(
os.path.join(
cfg["sat_image_path"],
str(this_year + 1),
"*h{:02d}v{:02d}*.h*".format(h, v),
)
)
)[:8]

inp_files = np.concatenate([inp_files1, inp_files2, inp_files3])
inp_files = np.array([i for i in inp_files if '.xml' not in i])
_, inds = np.unique([os.path.basename(f) for f in inp_files],
return_index=True)
inp_files = np.array([i for i in inp_files if ".xml" not in i])
_, inds = np.unique([os.path.basename(f) for f in inp_files], return_index=True)
inp_files = inp_files[inds]

vprm_inst = vprm(n_cpus=1, sites=[flux_tower_inst])

for c, i in enumerate(inp_files):
print(i)
if cfg['satellite'] == 'modis':
logger.info(i)
if cfg["satellite"] == "modis":
handler = modis(sat_image_path=i)
handler.load()
handler.crop(lonlat=(lon, lat), radius=4)
vprm_inst.add_sat_img(handler, b_nir='sur_refl_b02',
b_red='sur_refl_b01',
b_blue='sur_refl_b03',b_swir='sur_refl_b06',
which_evi='evi',
timestamp_key='sur_refl_day_of_year',
mask_bad_pixels=True,
mask_clouds=True,
mask_snow=True,
drop_bands=True)
vprm_inst.add_sat_img(
handler,
b_nir="sur_refl_b02",
b_red="sur_refl_b01",
b_blue="sur_refl_b03",
b_swir="sur_refl_b06",
which_evi="evi",
timestamp_key="sur_refl_day_of_year",
mask_bad_pixels=True,
mask_clouds=True,
mask_snow=True,
drop_bands=True,
)
else:
handler = VIIRS(sat_image_path=i)
handler.load()
handler.crop(lonlat=(lon, lat), radius=4)
vprm_inst.add_sat_img(handler, b_nir='SurfReflect_I2', b_red='SurfReflect_I1',
b_blue='no_blue_sensor', b_swir='SurfReflect_I3',
which_evi='evi2',
drop_bands=True)
vprm_inst.add_sat_img(
handler,
b_nir="SurfReflect_I2",
b_red="SurfReflect_I1",
b_blue="no_blue_sensor",
b_swir="SurfReflect_I3",
which_evi="evi2",
drop_bands=True,
)


# Uncomment if you want to do a spatial smearing with size=(3,3). Otherwise give size manually.

# Sort the satellite images by time and merge internally for easier computations
vprm_inst.sort_and_merge_by_timestamp()
vprm_inst.smearing(lonlats=[(lon, lat)],
keys=['evi', 'lswi'],
kernel=Gaussian2DKernel(1)) #(21, 21))
vprm_inst.smearing(
lonlats=[(lon, lat)], keys=["evi", "lswi"], kernel=Gaussian2DKernel(1)
) # (21, 21))

vprm_inst.reduce_along_lat_lon()

# Run lowess smoothing
vprm_inst.lowess(times=[datetime(this_year, 1, 1)+ timedelta(days=i) for i in np.arange(365.)],
frac=0.15, it=3,
keys=['evi', 'lswi'])
vprm_inst.lowess(
times=[datetime(this_year, 1, 1) + timedelta(days=i) for i in np.arange(365.0)],
frac=0.15,
it=3,
keys=["evi", "lswi"],
)

# Calculate necessary parameters for the vprm calculation
vprm_inst.calc_min_max_evi_lswi()
Expand All @@ -137,10 +177,12 @@ def all_files_exist(item):

# Loop over the year and get the vprm_variables
data_list = vprm_inst.data_for_fitting()
opath = '/home/b/b309233/software/VPRM_preprocessor/analysis_scripts/site_data_for_fit_modis/{}_{}.pickle'.format(this_year, s)
opath = "/home/b/b309233/software/VPRM_preprocessor/analysis_scripts/site_data_for_fit_modis/{}_{}.pickle".format(
this_year, s
)

if os.path.exists(opath):
os.remove(opath)
with open(opath, 'wb+') as ofile:

with open(opath, "wb+") as ofile:
pickle.dump(data_list, ofile)
Loading