willrayeo

#convert VV DN to Decibels
vv_dn = s1_cube.VV

vv_db = 10 * (xr.ufuncs.log10(vv_dn))
vv_db = vv_db.where(xr.ufuncs.isfinite(vv_db), 0)

# attribute long name and units to new variable
vv_db.attrs['long_name']='VV_dB'
vv_db.attrs['units']='decibels'

willrayeo

# assign all pixels equal or smaller than -20 a value of 1 and preserve the values of pixels 
step1 = s1_cube.VV_dB.where(s1_cube.VV_dB >= -20, 1)

# assign all other pixels a value of 0. 
flooded = step1.where(step1 == 1, 0)

flooded.attrs['long_name']='flooded'
flooded.attrs['units']='nounits'

s1_cube['flooded']= flooded

willrayeo

# Setup xcube
s2_cube_config = CubeConfig(dataset_name='S2L2A',
                         band_names=['B02', 'B03', 'B04', 'B08'],
                         bbox=tewkesbury_uk_bbox,
                         spatial_res=0.000089,
                         time_range=['2020-12-01', '2021-02-28'],
                         time_tolerance='30M')

# Open cube (on the fly)
s2_cube = open_cube(s2_cube_config, **sh_credentials)

willrayeo

timeseriesVV_dB = water_meadow.VV_dB.mean(axis=(1,2), skipna=True)
timeseriesVV_dB = timeseriesVV_dB.where(timeseriesVV_dB !=0).dropna("time")
timeseriesVV_dB

willrayeo

flood_sum = s1_cube.flooded.sum(dim="time")
flood_count = s1_cube.flooded.count(dim="time")
flood_average = flood_sum/flood_count

flood_sum.attrs['long_name']='flood_sum'
flood_sum.attrs['units']='nounits'

s1_cube['flood_sum']= flood_sum

flood_average.attrs['long_name']='flood_average'

willrayeo

# define NDWI in visualisation

ndwi = ((s2_cube.B03-s2_cube.B08)/(s2_cube.B03+s2_cube.B08))

ndwi.attrs['long_name']='NDWI'
ndwi.attrs['units']='unitless'

s2_cube['NDWI']= ndwi  

willrayeo

# define projection
flood_average_wgs84 = flood_average.rio.write_crs("epsg:4326")

# write to raster
flood_average_wgs84.rio.to_raster("flood_average.tif")

willrayeo

# define the timestep you wish to export
timestep = flood_threshold2_step2.sel(time='2020-12-28 10:00:00', method='nearest')

# define projection
timestep_wgs84 = timestep.rio.write_crs("epsg:4326")

# write to raster
timestep_wgs84.rio.to_raster("flood_extent_20201228.tif")

willrayeo

# Request the administrative boundaries as a bytes dataset.
france_admin_url = "https://www.data.gouv.fr/en/datasets/r/17062524-991f-4e13-9bf0-b410cc2216fd"
france_admin = requests.get(france_admin_url, stream=True)
france_admin_bytes = bytes(france_admin.content)

# Import administrative boundaries into a Geodataframe
with fiona.BytesCollection(france_admin_bytes) as f:
    crs = f.crs
    france_df = gpd.GeoDataFrame.from_features(f, crs=crs)
    

willrayeo

//VERSION=3
function setup() {
  return {
    input: ["B02", "B04", "B08", "B11"],
    output: [{id: "ndvi_t1", bands: 1, sampleType: "FLOAT32"},
             {id: "ndvi_t2", bands: 1, sampleType: "FLOAT32"},
             {id: "bsi_t2", bands: 1, sampleType: "FLOAT32"},
             {id: "harvested", bands: 1, sampleType: "UINT8"}],
    mosaicking: "ORBIT"
  };