Processing sar data

In this notebook, we will demonstrate the core functionality of the ReadSwesarr class. This class provides essential utilities for converting tif image to rasters and dataframe.

In [23]:

import sys
sys.path.append("../")
from swesarr_pytools.data_tools import ReadSwesarr, combine_swesarr_lidar
import pprint

import sys sys.path.append("../") from swesarr_pytools.data_tools import ReadSwesarr, combine_swesarr_lidar import pprint

Let's specify the directories we need to work it

In [12]:

# Just one band to test the downloading a singular band
x_band = "../resources/swesarr_data/GRMCT1_13801_19005_010_191106_225_XX_01/GRMCT1_13801_19005_010_191106_09225VH_XX_01.tif"

# Path to  v1 flight path
flight_folder = "../resources/swesarr_data/GRMCT1_13801_19005_010_191106_225_XX_01"

# Path to  v1 flight path
swesarr_fall_v1 = "../resources/swesarr_data/GRMCT1_13801_19005_010_191106_225_XX_01"
swesarr_winter_v1 = "../resources/swesarr_data/GRMCT2_13802_20006_012_200210_225_XX_01"

# Just one band to test the downloading a singular band x_band = "../resources/swesarr_data/GRMCT1_13801_19005_010_191106_225_XX_01/GRMCT1_13801_19005_010_191106_09225VH_XX_01.tif" # Path to v1 flight path flight_folder = "../resources/swesarr_data/GRMCT1_13801_19005_010_191106_225_XX_01" # Path to v1 flight path swesarr_fall_v1 = "../resources/swesarr_data/GRMCT1_13801_19005_010_191106_225_XX_01" swesarr_winter_v1 = "../resources/swesarr_data/GRMCT2_13802_20006_012_200210_225_XX_01"

We can also view a list of data files in a specified directory using the following command. These will be combined to form a 6D raster

In [6]:

example = ReadSwesarr(flight_path=flight_folder, version="v1")
example_list = example.get_data_files()
pprint.pprint(f"Folder List: {example_list}")

example = ReadSwesarr(flight_path=flight_folder, version="v1") example_list = example.get_data_files() pprint.pprint(f"Folder List: {example_list}")

("Folder List: ['GRMCT1_13801_19005_010_191106_13225VH_XX_01.tif', "
 "'GRMCT1_13801_19005_010_191106_17225VH_XX_01.tif', "
 "'GRMCT1_13801_19005_010_191106_09225VV_XX_01.tif', "
 "'GRMCT1_13801_19005_010_191106_17225VV_XX_01.tif', "
 "'GRMCT1_13801_19005_010_191106_13225VV_XX_01.tif', "
 "'GRMCT1_13801_19005_010_191106_09225VH_XX_01.tif']")

To retrieve the raster run the following . The resulting xarray dataset, can be manipulated using normal rioxarray/xarray operations, Since it is an xarray object itself.

In [8]:

example_raster = example.get_swesarr_raster()
example_raster

example_raster = example.get_swesarr_raster() example_raster

Out[8]:

<xarray.DataArray (band: 6, y: 3651, x: 2963)>
dask.array<concatenate, shape=(6, 3651, 2963), dtype=float32, chunksize=(1, 1200, 1200), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) <U4 '13VH' '17VH' '09VV' '17VV' '13VV' '09VH'
  * x            (x) float64 7.412e+05 7.412e+05 ... 7.471e+05 7.471e+05
  * y            (y) float64 4.328e+06 4.328e+06 4.328e+06 ... 4.32e+06 4.32e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    _FillValue:     nan
    scale_factor:   1.0
    add_offset:     0.0

xarray.DataArray

• band: 6
• y: 3651
• x: 2963

• dask.array<chunksize=(1, 1200, 1200), meta=np.ndarray>

  Array Chunk Bytes 247.60 MiB 5.49 MiB Shape (6, 3651, 2963) (1, 1200, 1200) Dask graph 72 chunks in 17 graph layers Data type float32 numpy.ndarray

• Coordinates: (4)
  • band
    (band)
    <U4
    '13VH' '17VH' ... '13VV' '09VH'

    array(['13VH', '17VH', '09VV', '17VV', '13VV', '09VH'], dtype='<U4')

  • x
    (x)
    float64
    7.412e+05 7.412e+05 ... 7.471e+05

    array([741213.109977, 741215.109977, 741217.109977, ..., 747133.109977,
           747135.109977, 747137.109977])

  • y
    (y)
    float64
    4.328e+06 4.328e+06 ... 4.32e+06

    array([4327758.006409, 4327756.006409, 4327754.006409, ..., 4320462.006409,
           4320460.006409, 4320458.006409])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 12N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32612"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 12N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -111.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 12N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32612"]]

    GeoTransform :

    741212.1099766805 2.0 0.0 4327759.006409172 0.0 -2.0

    array(0)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index(['13VH', '17VH', '09VV', '17VV', '13VV', '09VH'], dtype='object', name='band'))

  • x
    PandasIndex

    PandasIndex(Index([741213.1099766805, 741215.1099766805, 741217.1099766805,
           741219.1099766805, 741221.1099766805, 741223.1099766805,
           741225.1099766805, 741227.1099766805, 741229.1099766805,
           741231.1099766805,
           ...
           747119.1099766805, 747121.1099766805, 747123.1099766805,
           747125.1099766805, 747127.1099766805, 747129.1099766805,
           747131.1099766805, 747133.1099766805, 747135.1099766805,
           747137.1099766805],
          dtype='float64', name='x', length=2963))

  • y
    PandasIndex

    PandasIndex(Index([4327758.006409172, 4327756.006409172, 4327754.006409172,
           4327752.006409172, 4327750.006409172, 4327748.006409172,
           4327746.006409172, 4327744.006409172, 4327742.006409172,
           4327740.006409172,
           ...
           4320476.006409172, 4320474.006409172, 4320472.006409172,
           4320470.006409172, 4320468.006409172, 4320466.006409172,
           4320464.006409172, 4320462.006409172, 4320460.006409172,
           4320458.006409172],
          dtype='float64', name='y', length=3651))

• Attributes: (4)

  AREA_OR_POINT :

  Area

  _FillValue :

  nan

  scale_factor :

  1.0

  add_offset :

  0.0

To retrieve a dataframe of the xarray above run the following. The resulting dataframe, can be manipulated using normal pandas operations.

In [10]:

example_df = example.get_swesarr_df()
example_df.tail(5)

example_df = example.get_swesarr_df() example_df.tail(5)

Out[10]:

	y	x	F09VV	F09VH	F13VV	F13VH	F17VV	F17VH
10589389	4.320612e+06	746393.109977	-17.756666	-30.629366	-31.941320	-32.845257	-40.240311	-38.687332
10589390	4.320612e+06	746395.109977	-17.891167	-30.814734	-32.145603	-33.278454	-38.791313	-41.042126
10589391	4.320612e+06	746397.109977	-17.743111	-30.678053	-33.162926	-33.813126	-38.793060	-40.096455
10592352	4.320610e+06	746393.109977	-17.801363	-30.358782	-29.526028	-33.857082	-38.059433	-42.774803
10592353	4.320610e+06	746395.109977	-18.255501	-29.933887	-33.326679	-29.271967	-38.081070	-39.193199

To retrieve the a single band specify the band as shown below.

In [11]:

x_example = ReadSwesarr(flight_path=x_band, band='x')
x_example_raster = x_example.get_swesarr_raster()
x_example_raster

x_example = ReadSwesarr(flight_path=x_band, band='x') x_example_raster = x_example.get_swesarr_raster() x_example_raster

Out[11]:

<xarray.DataArray (band: 1, y: 3651, x: 2963)>
dask.array<open_rasterio-14a7a97d22cd3240f4cb28ae0ea87e94<this-array>, shape=(1, 3651, 2963), dtype=float32, chunksize=(1, 1200, 1200), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) <U4 '09VH'
  * x            (x) float64 7.412e+05 7.412e+05 ... 7.471e+05 7.471e+05
  * y            (y) float64 4.328e+06 4.328e+06 4.328e+06 ... 4.32e+06 4.32e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    _FillValue:     nan
    scale_factor:   1.0
    add_offset:     0.0

xarray.DataArray

• band: 1
• y: 3651
• x: 2963

• dask.array<chunksize=(1, 1200, 1200), meta=np.ndarray>

  Array Chunk Bytes 41.27 MiB 5.49 MiB Shape (1, 3651, 2963) (1, 1200, 1200) Dask graph 12 chunks in 2 graph layers Data type float32 numpy.ndarray

• Coordinates: (4)
  • band
    (band)
    <U4
    '09VH'

    array(['09VH'], dtype='<U4')

  • x
    (x)
    float64
    7.412e+05 7.412e+05 ... 7.471e+05

    array([741213.109977, 741215.109977, 741217.109977, ..., 747133.109977,
           747135.109977, 747137.109977])

  • y
    (y)
    float64
    4.328e+06 4.328e+06 ... 4.32e+06

    array([4327758.006409, 4327756.006409, 4327754.006409, ..., 4320462.006409,
           4320460.006409, 4320458.006409])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 12N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32612"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 12N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -111.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 12N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32612"]]

    GeoTransform :

    741212.1099766805 2.0 0.0 4327759.006409172 0.0 -2.0

    array(0)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index(['09VH'], dtype='object', name='band'))

  • x
    PandasIndex

    PandasIndex(Index([741213.1099766805, 741215.1099766805, 741217.1099766805,
           741219.1099766805, 741221.1099766805, 741223.1099766805,
           741225.1099766805, 741227.1099766805, 741229.1099766805,
           741231.1099766805,
           ...
           747119.1099766805, 747121.1099766805, 747123.1099766805,
           747125.1099766805, 747127.1099766805, 747129.1099766805,
           747131.1099766805, 747133.1099766805, 747135.1099766805,
           747137.1099766805],
          dtype='float64', name='x', length=2963))

  • y
    PandasIndex

    PandasIndex(Index([4327758.006409172, 4327756.006409172, 4327754.006409172,
           4327752.006409172, 4327750.006409172, 4327748.006409172,
           4327746.006409172, 4327744.006409172, 4327742.006409172,
           4327740.006409172,
           ...
           4320476.006409172, 4320474.006409172, 4320472.006409172,
           4320470.006409172, 4320468.006409172, 4320466.006409172,
           4320464.006409172, 4320462.006409172, 4320460.006409172,
           4320458.006409172],
          dtype='float64', name='y', length=3651))

• Attributes: (4)

  AREA_OR_POINT :

  Area

  _FillValue :

  nan

  scale_factor :

  1.0

  add_offset :

  0.0

Next we will combined a two swesarr flight, reprojected to same coordinates with lidar snow depth data. Note the data version is specified in this call since we are using version three data.

In [24]:

# path to lidar snow depth measurement for the area
lidar_path = "../resources/lidar_data/ASO_GrandMesa_2020Feb13_snowdepth_3m.tif"

# Path to  v3 flight path
swesarr_fall_v3 = "../resources/swesarr_data_v3/19004_011"
swesarr_winter_v3 = "../resources/swesarr_data_v3/20006_003"

# path to lidar snow depth measurement for the area lidar_path = "../resources/lidar_data/ASO_GrandMesa_2020Feb13_snowdepth_3m.tif" # Path to v3 flight path swesarr_fall_v3 = "../resources/swesarr_data_v3/19004_011" swesarr_winter_v3 = "../resources/swesarr_data_v3/20006_003"

In [25]:

combined_df = combine_swesarr_lidar(
    fall_flight_directory=swesarr_fall_v3,
    winter_flight_directory=swesarr_winter_v3,
    lidar_flight_path=lidar_path,
    drop_na=False,
    version="v3")

combined_df = combine_swesarr_lidar( fall_flight_directory=swesarr_fall_v3, winter_flight_directory=swesarr_winter_v3, lidar_flight_path=lidar_path, drop_na=False, version="v3")