In [1]:
import numpy as np
%matplotlib inline
In [2]:
from datetime import datetime
# That I can sqap areas
vietnam_extent = dict(
lat = (10.9348, 11.0190),
lon = (107.8164, 107.9168),
)
menindee_extent = dict(
lat = ( -32.420,-32.272),
lon = (142.2348, 142.407),
)
kenya_extent = dict(
lat = (2.4787, 2.7887),
lon = (35.7330, 37.4963)
)
extent = kenya_extent
A quick visualization of the area before it is loaded through the datacube
In [3]:
from utils.data_cube_utilities.dc_display_map import display_map
display_map(latitude = extent['lat'],
longitude = extent['lon'])
Out[3]:
In [4]:
import datacube
dc = datacube.Datacube()
products = dc.list_products()
products[products['name'].str.contains("alos")]
Out[4]:
In [5]:
menindee_extent["product"] = "s1_gamma0_menindee_lakes"
menindee_extent["platform"] = "SENTINEL_1"
kenya_extent["product"] = "alos_palsar_kenya"
kenya_extent["platform"] = "ALOS"
In [6]:
from utils.data_cube_utilities.data_access_api import DataAccessApi
dca = DataAccessApi()
dca.get_datacube_metadata(platform = extent["platform"],
product = extent["product"])
In [7]:
product_details = dc.list_products()[dc.list_products().name ==extent["product"]]
product_details
Out[7]:
In [8]:
import datacube
dc = datacube.Datacube()
In [9]:
dc.list_measurements().loc[extent['product']]
Out[9]:
In [10]:
print(extent)
dataset = dc.load(**extent, measurements=['hh', 'hv'])
In [11]:
dataset
Out[11]:
In [12]:
first_slice = dataset.isel(time = 1) # iselect selects by index, rather than value.
In [13]:
first_slice
Out[13]:
Plotting helper
In [14]:
def figure_ratio(ds, fixed_width = 20):
width = fixed_width
height = len(ds.latitude) * (fixed_width / len(ds.longitude))
return (width, height)
In [15]:
%matplotlib inline
first_slice.hh.plot(cmap = "Greys",figsize = figure_ratio(first_slice))
Out[15]:
In [16]:
first_slice.hv.plot(cmap = "Greys", figsize = figure_ratio(first_slice))
Out[16]:
In [17]:
import xarray as xr
def augment_dataset_with_ratio(ds: xr.Dataset, band_name_1: str, band_name_2 :str) -> xr.Dataset:
a_per_b = (ds[band_name_1]/ds[band_name_2])
a_per_b = a_per_b.to_dataset(name = '{b1}_per_{b2}'.format(b1 = band_name_1, b2 = band_name_2)) # turn xarray.dataarray into xarray dataset and name dataset variable.
return ds.merge(a_per_b)
Apply function and display new xarray.Dataset
In [18]:
first_slice = augment_dataset_with_ratio(first_slice, "hv", "hh")
first_slice = augment_dataset_with_ratio(first_slice, "hh", "hv")
print(first_slice)
In [ ]:
np.log10(first_slice.hv_per_hh).plot(cmap = "Greys", figsize = figure_ratio(first_slice))
Out[ ]:
The function below defines a basic normaization and plotting function for xarrays. I would not recommend re-using this unless you understand your data in great enough detail to determine that this suites your needs. This is a first pass at including a plotting utility in our S1 analysis workflow. The canonical approach or desired results would be something along the lines of following established processing methods used in something like ESA's Sentinel-1 Toolbox or from the Sentinel Callibration Guide.
In [ ]:
import numpy as np
import xarray as xr
def build_rgb_from_ds(_dataset: xr.Dataset,
r:str = None,
g:str = None,
b:str = None,
logarithmic_scale_enabled = False):
norm = _dataset.copy()
if logarithmic_scale_enabled == True:
norm[r] = np.log10(norm[r])
norm[g] = np.log10(norm[g])
norm[b] = np.log10(norm[b])
norm = (255 * _dataset/_dataset.max()).astype(np.uint16)
_r = norm[r].values #.astype(np.float32)
_g = norm[g].values #.astype(np.float32)
_b = norm[b].values #_per_vh.astype(np.float32)
_rgb = np.dstack([_r,_g,_b])
return _rgb
In [ ]:
rgb = build_rgb_from_ds(first_slice,
r = "hv",
g = "hh",
b = "hv_per_hh")
In [ ]:
import matplotlib.pyplot as plt
plt.figure( figsize = figure_ratio(first_slice))
plt.imshow( rgb )
Out[ ]:
In [ ]:
from utils.data_cube_utilities.dc_chunker import create_geographic_chunks, combine_geographic_chunks
geographic_chunks = create_geographic_chunks(latitude=extent["lat"], longitude=extent["lon"], geographic_chunk_size=.05)
In [ ]:
for x in geographic_chunks:
print(x)
In [ ]:
from utils.data_cube_utilities.dc_mosaic import create_median_mosaic
from utils.data_cube_utilities.dc_sar_utils import dn_to_db
import warnings
warnings.filterwarnings("ignore")
import numpy as np
measurements = ['hv', 'hh']
product_chunks = []
# This part is unpythonic but explicit
for index, chunk in enumerate(geographic_chunks):
data = dca.get_dataset_by_extent(extent["product"],
longitude=chunk['longitude'], latitude=chunk['latitude'],
measurements=measurements)
if 'hv' in data:
product_chunks.append(create_median_mosaic(data, dtype="float32", no_data=0))
final_mosaic = combine_geographic_chunks(product_chunks)
In [ ]:
final_mosaic
Out[ ]:
In [ ]:
final_mosaic = augment_dataset_with_ratio(final_mosaic, "hv", "hh")
final_mosaic = augment_dataset_with_ratio(final_mosaic, "hh", "hv")
In [ ]:
final_mosaic
Out[ ]:
In [ ]:
rgb = build_rgb_from_ds(final_mosaic,
r = "hv",
g = "hh",
b = "hh_per_hv")
In [ ]:
import matplotlib.pyplot as plt
plt.figure(figsize = figure_ratio(final_mosaic))
plt.imshow( rgb )
Out[ ]:
In [ ]:
dataset = augment_dataset_with_ratio(dataset, "hv", "hh") ## Compute HV/HH for entire dataset rather than just a slice(see section above)
In [ ]:
display_map(latitude = extent['lat'], longitude = extent['lon'])
In [ ]:
# Lat and Lon coordinates extracted from the map above
pixel_lat = -32.3626 #Menindee
pixel_lon = 142.2764 #Menindee
pixel_lat = 11.1306
pixel_lon = 107.6052
Select a pixel from our xarray.Dataset using nearest neighbor.
In [ ]:
pixel = dataset.sel(latitude = pixel_lat,
longitude = pixel_lon,
method = 'nearest') # nearest neighbor selection
In [ ]:
pixel
In [ ]:
pixel.hh.plot.hist()
In [ ]:
pixel.hv.plot.hist()
In [ ]:
## X-axis time, Y-axis values of vh
plt.figure(figsize = (20,5))
plt.scatter(pixel.time.values, pixel.hh.values)
In [ ]:
## X-axis time, Y-axis values of vh
plt.figure(figsize = (20,5))
plt.scatter(pixel.time.values, np.log10(pixel.hh.values))
In [ ]:
## X-axis time, Y-axis values of vv
plt.figure(figsize = (20,5))
plt.scatter(pixel.time.values, pixel.hv.values)
In [ ]:
## X-axis time, Y-axis values of vv
plt.figure(figsize = (20,5))
plt.scatter(pixel.time.values, np.log10(pixel.hv.values))
In [ ]:
## X-axis time, Y-axis values of vv
plt.figure(figsize = (20,5))
plt.scatter(pixel.time.values, pixel.hv_per_hh.values)
In [ ]:
## X-axis time, Y-axis values of vv
plt.figure(figsize = (20,5))
plt.scatter(pixel.time.values,
(np.log10(pixel.hv) / np.log10(pixel.hh)).values
)
In [ ]:
def ds_to_timeseries(ds, band, log = False):
da = ds[band]
da = np.log10(da) if log == True else da
da_by_time_slices = [da.isel(time = i) for i in range(len(da.time))]
# return [x.values.flatten() for x in da_by_time_slices]
return [x.values[np.isfinite(x.values)].flatten() for x in da_by_time_slices]
In [ ]:
plt.figure(figsize = (20,5))
plt.boxplot(ds_to_timeseries(dataset, "hv"), 0, "");
In [ ]:
plt.figure(figsize = (20,5))
plt.boxplot(ds_to_timeseries(dataset, "hh"), 0, "");
In [ ]:
plt.figure(figsize = (20,5))
plt.boxplot(ds_to_timeseries(dataset, "hv_per_hh"), 0, "");
In [ ]:
import matplotlib
def plot_threshold(ds, band_name, bottom = None , top = None, log_scaled = False, cmap_name = 'Greys'):
# Threshold is applied to original data, not log scaled data(if you haven't scaled already)
_range = "Full {} range: {}-{}".format(band_name, ds[band_name].min().values,ds[band_name].max().values )
def figure_ratio(ds, fixed_width = 20):
width = fixed_width
height = len(ds.latitude) * (fixed_width / len(ds.longitude))
return (width, height)
selection = ds[band_name]
my_cmap = matplotlib.cm.get_cmap(cmap_name)
my_cmap.set_over('r')
my_cmap.set_under('b')
plt.figure(figsize = figure_ratio(ds))
selection = np.log10(selection) if log_scaled == True else selection
bottom = np.log10(bottom) if log_scaled == True and bottom is not None else bottom
top = np.log10(top) if log_scaled == True and top is not None else top
selection.plot(cmap = my_cmap, vmax =top, vmin = bottom)
plt.figtext(0.7,0,_range, horizontalalignment = "center")
In [ ]:
plot_threshold(first_slice, "hh", log_scaled = False)
In [ ]:
plot_threshold(first_slice, "hv",
bottom = 437,
top = 3500)
In [ ]:
plot_threshold(first_slice, "hh")
In [ ]:
plot_threshold(first_slice,
"hh",
bottom = 1000,
top = 10000
)