# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Crop-type image library
#
# Using example locations of crop types in different regions of Africa, create example images and NDVI timeseries to assist in the training and validation data collections
#
# Image library is here: https://docs.google.com/document/d/1lLOYIM12_E09kC6l0ezUzBngYFTC7DDEIftAZzmvn2M/edit?usp=sharing

import geopandas as gpd
import datacube
import numpy as np
import sys
import xarray as xr
from datacube.utils import geometry
import matplotlib.pyplot as plt 
sys.path.append('../Scripts')
from deafrica_datahandling import load_ard
from deafrica_bandindices import calculate_indices
from deafrica_plotting import rgb, map_shapefile
from deafrica_spatialtools import xr_rasterize
from deafrica_classificationtools import HiddenPrints

dc = datacube.Datacube(app='image_lib')

# ## Analysis Parameters

shp_file= 'data/crop_type_library/madagascar/Itasy_rice.shp'
product = 's2_l2a'#'usgs_ls8c_level2_2'# 'ls8_usgs_sr_scene'
time_ndvi = ('2018-06', '2020-06')
time_rgb = ('2019')

# ## Open file and check extents

gdf = gpd.read_file(shp_file).to_crs('epsg:4326')
gdf.head()

print(gdf['layer'].unique())

crop='sweet potatoes'

gdf = gdf[gdf['Comment'] == crop]
gdf

gdf = gdf.sample(n=3).reset_index()

for index,row in gdf.iterrows():
    gdf.loc[index, 'Crop'] = crop+str(index)

map_shapefile(gdf, attribute='Crop', hover_col='Crop')

# ## Extract RGB and NDVI time-series

# +
# gdf = gdf.to_crs('EPSG:6933')
# gdf['geometry'] = gdf['geometry'].buffer(20).envelope
# gdf = gdf.to_crs('EPSG:4326')

# +
ndvi_dict = {}
rgb_dict = {}
for index, row in gdf.iterrows():
    name=gdf.iloc[[index]].Crop.values[0]
#     dates = gdf.iloc[[index]]['Time'].values[0]
#     print(name + ": "+ str(dates))
    geom = geometry.Geometry(
            gdf.geometry.values[index].__geo_interface__, geometry.CRS(
                    'epsg:4326'))
    query_ndvi = {
        "geopolygon": geom,
        'time': time_ndvi,
        'measurements': ['red',
                         'nir'],
        'resolution': (-10, 10),
        'group_by': 'solar_day',
        'output_crs':'EPSG:6933'
    }

    # load cloud-masked fractional cover using load_ard
    with HiddenPrints():
        ds = load_ard(dc=dc,
                  **query_ndvi,
                  products=[product], 
                  min_gooddata=0.50
                 )
    
    with HiddenPrints():
        mask = xr_rasterize(gdf.iloc[[index]], ds)
    
#     ds = ds.where(mask)
    ds = calculate_indices(ds, 'NDVI', collection='s2')
    ds_index = ds.mean(['x','y'])
    ds_index = ds_index.interpolate_na(dim="time", method="linear")
    ds_index = ds_index.resample(time="1M").median()
    ndvi_dict[name] = ds_index
    
#     query_rgb = {
#         "geopolygon": geom,
#         'time': time_rgb,
#         'measurements': ['green',
#                          'blue',
#                          'red'],
#         'resolution': (-10, 10),
#         'group_by': 'solar_day',
#         'output_crs':'EPSG:6933'
#     }
    
#     # load cloud-masked fractional cover using load_ard
#     with HiddenPrints():
#         ds = load_ard(dc=dc,
#                   **query_rgb,
#                   products=[product],
#                   min_gooddata=1.0
#                  )
    
#     with HiddenPrints():
#         mask = xr_rasterize(gdf.iloc[[index]], ds)
    
#     ds = ds.where(mask)
    
#     rgb_dict[name]=ds
# -

key = crop+'1'

ndvi_dict[key].NDVI.plot.line('b--', figsize=(11,5))
plt.ylim(0,1.0)
plt.title(key[0:-1] +" "+ product)

# +
# rgb(rgb_dict[key], index=list(range(0,len(rgb_dict[key].time))))
# -



StarcoderdataJupyterScriptsDedupFiltered