Saroumanedraft example

try:
    client.close()
    cluster.close()
except:
    pass
import bokeh.io
bokeh.io.output_notebook()
import holoviews as hv
hv.extension('bokeh')
import datashader as ds
from holoviews.operation.datashader import datashade,rasterize
from IPython.display import display, HTML

Loading BokehJS ...

import socket,os
hostname = socket.gethostname()
from dask.distributed import Client, LocalCluster, progress
cluster = LocalCluster(n_workers=4,
                       threads_per_worker=4,
                       scheduler_port=0,
                       host=hostname,
                       dashboard_address='%s:0' % hostname,
                       memory_limit='8GB',
                       local_directory=os.path.join("/tmp",os.getenv("USER"),"dask")
                      )
client = Client(cluster)
client

Client Scheduler: tcp://127.0.1.1:42926 Dashboard: http://127.0.1.1:39919/status

Cluster Workers: 4 Cores: 32 Memory: 32.00 GB

import numpy as np
import saroumanedraft as sr
import logging
sr.logger.setLevel(logging.DEBUG)

chunks={'pol' : 1 , 'xtrack' : 500 , 'atrack' : 200  }
sar_ds = sr.open_dataset('/home/datawork-cersat-public/cache/project/mpc-sentinel1/data/esa/sentinel-1a/L1/IW/S1A_IW_GRDH_1S/2017/250/S1A_IW_GRDH_1SDV_20170907T103020_20170907T103045_018268_01EB76_992F.SAFE',chunks=chunks)
sar_ds

DEBUG:SentinelReader:S1 dataset opened/merged in 1.2 s

xarray.Dataset

• Dimensions:
  • atrack: 16778
  • pol: 2
  • xtrack: 25187
• Coordinates: (3)
  • pol
    (pol)
    object
    'VV' 'VH'

    array(['VV', 'VH'], dtype=object)

  • xtrack
    (xtrack)
    int64
    0 1 2 3 ... 25183 25184 25185 25186

    array([    0,     1,     2, ..., 25184, 25185, 25186])

  • atrack
    (atrack)
    int64
    0 1 2 3 ... 16774 16775 16776 16777

    array([    0,     1,     2, ..., 16775, 16776, 16777])

• Data variables: (4)
  • digital_number
    (pol, xtrack, atrack)
    uint16
    dask.array<chunksize=(1, 500, 200), meta=np.ndarray>

    
    
    
    
      
    
        
     
     Array 
     Chunk 
      
      
    
        
     Bytes 
     1.69 GB 
     
     200.00 kB 
        
     Shape 
     (2, 25187, 16778) 
     
     (1, 500, 200) 
        
     Count 
     42842 Tasks 
     8568 Chunks 
        
     Type 
     uint16 
     numpy.ndarray 
      
    
    
    
    
    
    
    

  • longitude
    (xtrack, atrack)
    float64
    dask.array<chunksize=(500, 200), meta=np.ndarray>

    
    
    
    
      
    
        
     
     Array 
     Chunk 
      
      
    
        
     Bytes 
     3.38 GB 
     
     800.00 kB 
        
     Shape 
     (25187, 16778) 
     
     (500, 200) 
        
     Count 
     30395 Tasks 
     4284 Chunks 
        
     Type 
     float64 
     numpy.ndarray 
      
    
    
    
    
    
    
    

  • latitude
    (xtrack, atrack)
    float64
    dask.array<chunksize=(500, 200), meta=np.ndarray>

    
    
    
    
      
    
        
     
     Array 
     Chunk 
      
      
    
        
     Bytes 
     3.38 GB 
     
     800.00 kB 
        
     Shape 
     (25187, 16778) 
     
     (500, 200) 
        
     Count 
     30395 Tasks 
     4284 Chunks 
        
     Type 
     float64 
     numpy.ndarray 
      
    
    
    
    
    
    
    

  • sigma0
    (pol, xtrack, atrack)
    float64
    dask.array<chunksize=(1, 500, 200), meta=np.ndarray>

    
    
    
    
      
    
        
     
     Array 
     Chunk 
      
      
    
        
     Bytes 
     6.76 GB 
     
     800.00 kB 
        
     Shape 
     (2, 25187, 16778) 
     
     (1, 500, 200) 
        
     Count 
     85682 Tasks 
     8568 Chunks 
        
     Type 
     float64 
     numpy.ndarray 
      
    
    
    
    
    
    
    

• Attributes: (9)

  product_class :

  S

  product_type :

  GRD

  mission :

  SENTINEL-1

  satellite :

  A

  instrument_mode :

  IW

  instrument_swath :

  IW

  start_date :

  2017-09-07 10:30:20.936409

  stop_date :

  2017-09-07 10:30:45.935264

  coordinates :

  [[19.215193, -68.158363], [19.640442, -70.514343], [21.147583, -70.221626], [20.725643, -67.842209]]

# take a slice : still no evaluation
#sliced_sar = sar_ds.sel(xtrack=slice(0,500),atrack=slice(150,300),pol='VV')
#sliced_sar = sar_ds.sel(xtrack=slice(0,500),atrack=slice(150,300),pol='VH')
#sliced_sar = sar_ds.sel(xtrack=slice(0,25186,100),atrack=slice(0,16777,100),pol='VH')
# eye
sliced_sar = sar_ds.sel(xtrack=slice(5000,9999,1),atrack=slice(7000,11999,1),pol='VH')
# slice with earth at high res
#sliced_sar = sar_ds.sel(xtrack=slice(24000,25186,1),atrack=slice(14500,15500,1),pol='VH')
sliced_sar

xarray.Dataset

• Dimensions:
  • atrack: 5000
  • xtrack: 5000
• Coordinates: (3)
  • pol
    ()
    <U2
    'VH'

    array('VH', dtype='<U2')

  • xtrack
    (xtrack)
    int64
    5000 5001 5002 ... 9997 9998 9999

    array([5000, 5001, 5002, ..., 9997, 9998, 9999])

  • atrack
    (atrack)
    int64
    7000 7001 7002 ... 11998 11999

    array([ 7000,  7001,  7002, ..., 11997, 11998, 11999])

• Data variables: (4)
  • digital_number
    (xtrack, atrack)
    uint16
    dask.array<chunksize=(500, 200), meta=np.ndarray>

    
    
    
    
      
    
        
     
     Array 
     Chunk 
      
      
    
        
     Bytes 
     50.00 MB 
     
     200.00 kB 
        
     Shape 
     (5000, 5000) 
     
     (500, 200) 
        
     Count 
     43092 Tasks 
     250 Chunks 
        
     Type 
     uint16 
     numpy.ndarray 
      
    
    
    
    
    
    
    

  • longitude
    (xtrack, atrack)
    float64
    dask.array<chunksize=(500, 200), meta=np.ndarray>

    
    
    
    
      
    
        
     
     Array 
     Chunk 
      
      
    
        
     Bytes 
     200.00 MB 
     
     800.00 kB 
        
     Shape 
     (5000, 5000) 
     
     (500, 200) 
        
     Count 
     30645 Tasks 
     250 Chunks 
        
     Type 
     float64 
     numpy.ndarray 
      
    
    
    
    
    
    
    

  • latitude
    (xtrack, atrack)
    float64
    dask.array<chunksize=(500, 200), meta=np.ndarray>

    
    
    
    
      
    
        
     
     Array 
     Chunk 
      
      
    
        
     Bytes 
     200.00 MB 
     
     800.00 kB 
        
     Shape 
     (5000, 5000) 
     
     (500, 200) 
        
     Count 
     30645 Tasks 
     250 Chunks 
        
     Type 
     float64 
     numpy.ndarray 
      
    
    
    
    
    
    
    

  • sigma0
    (xtrack, atrack)
    float64
    dask.array<chunksize=(500, 200), meta=np.ndarray>

    
    
    
    
      
    
        
     
     Array 
     Chunk 
      
      
    
        
     Bytes 
     200.00 MB 
     
     800.00 kB 
        
     Shape 
     (5000, 5000) 
     
     (500, 200) 
        
     Count 
     85932 Tasks 
     250 Chunks 
        
     Type 
     float64 
     numpy.ndarray 
      
    
    
    
    
    
    
    

• Attributes: (9)

  product_class :

  S

  product_type :

  GRD

  mission :

  SENTINEL-1

  satellite :

  A

  instrument_mode :

  IW

  instrument_swath :

  IW

  start_date :

  2017-09-07 10:30:20.936409

  stop_date :

  2017-09-07 10:30:45.935264

  coordinates :

  [[19.215193, -68.158363], [19.640442, -70.514343], [21.147583, -70.221626], [20.725643, -67.842209]]

# dn at high res
future = client.compute(sliced_sar.digital_number)
progress(future)

dn_high = future.result()
display(dn_high)
rasterize(hv.Image(dn_high).opts(cmap='gray'))

xarray.DataArray

'digital_number'

• xtrack: 5000
• atrack: 5000

• 98 112 98 92 99 84 62 52 62 92 ... 104 129 129 87 65 78 77 64 52 63

  array([[ 98, 112,  98, ...,  89,  92, 111],
         [ 81,  69,  52, ...,  92,  84, 108],
         [ 61,  55,  60, ...,  69,  64,  92],
         ...,
         [ 57,  62,  66, ...,  86,  59,  63],
         [ 40,  59,  61, ...,  72,  49,  73],
         [ 51,  63,  46, ...,  64,  52,  63]], dtype=uint16)

• Coordinates: (3)
  • pol
    ()
    <U2
    'VH'

    array('VH', dtype='<U2')

  • xtrack
    (xtrack)
    int64
    5000 5001 5002 ... 9997 9998 9999

    array([5000, 5001, 5002, ..., 9997, 9998, 9999])

  • atrack
    (atrack)
    int64
    7000 7001 7002 ... 11998 11999

    array([ 7000,  7001,  7002, ..., 11997, 11998, 11999])

• Attributes: (0)

# sigma0 at high res
future = client.compute(sliced_sar.sigma0)
progress(future)

sigma0_high = future.result()
display(sigma0_high)
rasterize(hv.Image(sigma0_high).opts(cmap='gray',colorbar=True,tools=['hover']))

xarray.DataArray

'sigma0'

• xtrack: 5000
• atrack: 5000

• 0.02407 0.03143 0.02407 0.02121 ... 0.01605 0.01109 0.007322 0.01075

  array([[0.02406504, 0.0314319 , 0.02406504, ..., 0.0198479 , 0.02120851,
          0.03087312],
         [0.01644038, 0.01192999, 0.00677561, ..., 0.02120887, 0.01768074,
          0.02922735],
         [0.00932414, 0.00758009, 0.00902094, ..., 0.01193019, 0.01026382,
          0.02120923],
         ...,
         [0.00879739, 0.01040849, 0.01179484, ..., 0.02002632, 0.00942558,
          0.01074695],
         [0.00433242, 0.00942572, 0.01007558, ..., 0.01403703, 0.00650133,
          0.01442966],
         [0.00704299, 0.01074726, 0.0057297 , ..., 0.01109115, 0.00732189,
          0.01074726]])

• Coordinates: (3)
  • pol
    ()
    <U2
    'VH'

    array('VH', dtype='<U2')

  • xtrack
    (xtrack)
    int64
    5000 5001 5002 ... 9997 9998 9999

    array([5000, 5001, 5002, ..., 9997, 9998, 9999])

  • atrack
    (atrack)
    int64
    7000 7001 7002 ... 11998 11999

    array([ 7000,  7001,  7002, ..., 11997, 11998, 11999])

• Attributes: (0)

# do some science : move a box of shape (20,20) across sigma0 array, and compute mean by box
def science_function(sigma0):
    # here, sigma0 is a simple numpy array
    return np.mean(sigma0),np.std(sigma0)  # so the output size will be spatialy reduced, but expanding axis in a new dim

# shape_out is optionnal. It's the shape of returned values from science_function, with the count of returned 
# values as the first dim. setting shape_out=None will compute this value for the user, but for better optimisation,
# the user should hard code the value.
sigma0_low,sigma0_std = sliced_sar.sarlib.map(science_function,sliced_sar.sigma0,shape_in=(20,10),shape_out=(2, 1, 1))
sigma0_low 

Array Chunk Bytes 1000.00 kB 4.00 kB Shape (250, 500) (25, 20) Count 87432 Tasks 250 Chunks Type float64 numpy.ndarray

# compute science
future = client.compute((sigma0_low,sigma0_std))
progress(future)

sigma0_low = future[0].result()
sigma0_std = future[1].result()

display(sigma0_low.shape)
rasterize(hv.Image(sigma0_low).opts(cmap='gray',colorbar=True,tools=['hover'])) + rasterize(hv.Image(sigma0_std).opts(cmap='gray',colorbar=True,tools=['hover']))

(250, 500)

# sigma0 debug
#from importlib import reload
#reload(sr)
#sr.logger.setLevel(logging.DEBUG)
#import sys
#sys.path.append('/home3/homedir7/perso/oarcher/bin/')
#import pydev
lut = sar_ds.sarlib.metadata._pol[sliced_sar.pol.item()]['sigma0_lut']
xtrack_val = np.array(dn_high.xtrack).copy()
atrack_val = np.array(dn_high.atrack).copy()
sigma0 = sr._get_sigma0_from_lut(
    np.asarray(dn_high).copy(),
    xtrack_val,
    atrack_val,
    np.asarray(lut).copy(),
    np.asarray(lut.xtrack).copy(),
    np.asarray(lut.atrack).copy(),
    block_info=((xtrack_val[0],xtrack_val[-1]),(atrack_val[0],atrack_val[-1]))
)
#rasterize(hv.Image(np.asarray(sigma0)).redim.range(z=(0, 1e-2))) + rasterize(hv.Image(np.asarray(dn_high)).opts(hv.opts.Image(tools=['hover'],colorbar=True)))
#[ sigma0 , np.array(dn_val) ]
#(np.asarray(dn_val),np.asarray(sliced_sar.xtrack),np.asarray(sliced_sar.atrack.data),np.asarray(lut),np.asarray(lut.xtrack),np.asarray(lut.atrack))
#hv.Image(np.asarray(sigma0))

DEBUG:SentinelReader:_get_sigma0_from_lut not called from map_blocks.
DEBUG:SentinelReader:block position : ((5000, 10000), (7000, 12000))

rasterize(hv.Image(np.asarray(sigma0)).opts(cmap='gray',colorbar=True,tools=['hover']))

# manual interp from lut

import xarray as xr
import time

t1 = time.time()

from scipy.interpolate import RectBivariateSpline
xx , yy = np.mgrid[xtrack_val[0]:xtrack_val[-1]+1,atrack_val[0]:atrack_val[-1]+1]
interp_lut = RectBivariateSpline(lut.xtrack,lut.atrack,lut,kx=1,ky=1).ev(xx,yy) 
interp_lut = xr.DataArray(interp_lut,dims=['xtrack','atrack'], coords = {'xtrack' : dn_high.xtrack.values , 'atrack' : dn_high.atrack.values})

#xarray interp is slowest than RectBivariateSpline
#interp_lut = lut.interp({'xtrack' : dn_high.xtrack , 'atrack' : dn_high.atrack},method='linear')

sigma0_1 = np.abs(dn_high)**2 / (interp_lut**2)
diff=np.max(np.abs(sigma0_1-sigma0))

print("computed in %2.1f secs. max diff : %f" % (time.time()-t1 , diff))

#(hv.HeatMap(lut)  +  hv.HeatMap(interp_lut)).opts(hv.opts.HeatMap(axiswise=True,tools=['hover'],colorbar=True)).redim.range(z=(500, 700))
#hv.Image(sigma0)

computed in 56.3 secs. max diff : 0.000000