Download and convert to Geopackage

This downloads SWOT Pixel Cloud products from hydroweb.next (API-Key necessary) based on a region and a period of interest. Then is extracts information contained in the area of interest for your study, stores everything in a Geopackage Database for future use. Geopackage is a convenient data storage format, based on SQL, and is compatible with QGIS.

Setting the region and period of interest

Using a geopackage layer, preliminary created with, e.g. QGIS, to limit data download and database

from pixcdust.downloaders.hydroweb_next import PixCDownloader
import geopandas as gpd
from datetime import datetime

# reading the area of interest
gdf_geom = gpd.read_file("../data/aoi.gpkg")

# Limiting time period
dates = (
    datetime(2023,4,6),
    datetime(2023,4,8),
)

Download

This will unfortunately lead to downloading many big files (that will be removed later). This is the only way right now, but the hydroweb.next team is working on improving that.

pixcdownloader = PixCDownloader(
    gdf_geom,
    dates,
    verbose=1,
    path_download='/tmp/pixc',
    )
pixcdownloader.search_download()

Extraction

Now we have all necessary files, let us extract key variables within area of interest in a geopackage database. This geopackage format is quite efficient (though not the most efficient), and may easily be visualized in, e.g., QGIS We are using the same geodataframe to limit the data to the area of interest

from pixcdust.converters.gpkg import Nc2GpkgConverter
from glob import glob

# You can specify conditions on variables to filter data
conditions= {"sig0":{'operator': "gt", 'threshold': 20},  # sig0 > 20
             "classification":{'operator': "ge", 'threshold': 3},  # classification >= 3
            }

pixc = Nc2GpkgConverter(
            path_in = glob(pixcdownloader.path_download+'/*/*nc'),
            variables=['height', 'sig0', 'classification'],
            area_of_interest=gdf_geom,
            conditions=conditions
        )
pixc.database_from_nc(path_out="/tmp/pixc_gpkg.gpkg")

100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 32.53it/s]

skipping layer 20230407_483_16_78L                             (already in geopackage /tmp/pixc_gpkg.gpkg)
skipping layer 20230406_482_16_78L                             (already in geopackage /tmp/pixc_gpkg.gpkg)

database has been succesfully created, we can remove the raw files

# import shutil
# shutil.rmtree('/tmp/pixc')

Read the database

Previous steps are not necessary

Now we can open this database in a GeoDataFrame, load it in, e.g., QGIS, etc.

from pixcdust.readers.gpkg import GpkgReader

# nb: you may specify
pixc_read = GpkgReader(
    "/tmp/pixc_gpkg.gpkg"
)
pixc_read.read()
pixc_read.data

100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 10.67it/s]

<xarray.Dataset> Size: 4MB
Dimensions:         (index: 49160)
Coordinates:
  * index           (index) int64 393kB 0 1 2 3 4 ... 49156 49157 49158 49159
Data variables:
    height          (index) float64 393kB 324.7 284.0 274.1 ... 165.7 170.7
    sig0            (index) float64 393kB 42.92 20.86 39.61 ... 24.67 22.21 28.0
    classification  (index) float64 393kB 6.0 3.0 3.0 3.0 ... 3.0 3.0 3.0 3.0
    geoid           (index) float64 393kB 49.44 49.44 49.42 ... 49.32 49.32
    latitude        (index) float64 393kB 43.52 43.52 43.54 ... 43.68 43.68
    longitude       (index) float64 393kB 1.462 1.459 1.459 ... 1.375 1.379
    wse             (index) float64 393kB 275.2 234.5 224.7 ... 116.4 121.3
    geometry        (index) geometry 393kB <class 'xarray.core.extension_arra...

xarray.Dataset

• Dimensions:
  • index: 49160
• Coordinates: (1)
  • index
    (index)
    int64
    0 1 2 3 ... 49156 49157 49158 49159

    array([    0,     1,     2, ..., 49157, 49158, 49159])

• Data variables: (8)
  • height
    (index)
    float64
    324.7 284.0 274.1 ... 165.7 170.7

    array([324.65289307, 283.96170044, 274.08056641, ..., 165.60694885,
           165.71679688, 170.66943359])

  • sig0
    (index)
    float64
    42.92 20.86 39.61 ... 22.21 28.0

    array([42.91566849, 20.8577652 , 39.60572052, ..., 24.66830826,
           22.21314812, 28.00420761])

  • classification
    (index)
    float64
    6.0 3.0 3.0 3.0 ... 3.0 3.0 3.0 3.0

    array([6., 3., 3., ..., 3., 3., 3.])

  • geoid
    (index)
    float64
    49.44 49.44 49.42 ... 49.32 49.32

    array([49.44259262, 49.43933487, 49.42419052, ..., 49.31718826,
           49.31721878, 49.31978607])

  • latitude
    (index)
    float64
    43.52 43.52 43.54 ... 43.68 43.68

    array([43.51782585, 43.52280247, 43.53983836, ..., 43.68376816,
           43.68379813, 43.68208302])

  • longitude
    (index)
    float64
    1.462 1.459 1.459 ... 1.375 1.379

    array([1.46168274, 1.45885209, 1.45871443, ..., 1.37530015, 1.37546673,
           1.37948554])

  • wse
    (index)
    float64
    275.2 234.5 224.7 ... 116.4 121.3

    array([275.21029663, 234.52236938, 224.65637207, ..., 116.2897644 ,
           116.39958191, 121.34964752])

  • geometry
    (index)
    geometry
    <class 'xarray.core.extension_ar...

    <class 'xarray.core.extension_array.PandasExtensionArray'>(array=<GeometryArray>
    [<POINT (1.462 43.518)>, <POINT (1.459 43.523)>,  <POINT (1.459 43.54)>,
     <POINT (1.463 43.543)>, <POINT (1.459 43.524)>, <POINT (1.458 43.534)>,
     <POINT (1.458 43.539)>, <POINT (1.458 43.539)>,  <POINT (1.458 43.54)>,
      <POINT (1.458 43.54)>,
     ...
     <POINT (1.371 43.679)>, <POINT (1.371 43.679)>, <POINT (1.372 43.679)>,
     <POINT (1.372 43.679)>,  <POINT (1.372 43.68)>, <POINT (1.372 43.682)>,
     <POINT (1.372 43.679)>, <POINT (1.375 43.684)>, <POINT (1.375 43.684)>,
     <POINT (1.379 43.682)>]
    Length: 49160, dtype: geometry)

• Indexes: (1)
  • index
    PandasIndex

    PandasIndex(RangeIndex(start=0, stop=49160, step=1, name='index'))

• Attributes: (0)

Rasterization into H3 grid

from pixcdust.converters.gpkg import GpkgDGGSProjecter
h3_grid = GpkgDGGSProjecter("/tmp/pixc_gpkg.gpkg", 10, path_out = '../data/h3_gpkg.gpkg', healpix=False) # True for healpix projection
h3_grid.compute_layers()

Layers: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:01<00:00,  1.42it/s]

pixc_read = GpkgReader(
    '../data/h3_gpkg.gpkg'
)
pixc_read.read()
pixc_read.data

100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 98.45it/s]

<xarray.Dataset> Size: 725kB
Dimensions:         (index: 11333)
Coordinates:
  * index           (index) int64 91kB 0 1 2 3 4 ... 11329 11330 11331 11332
Data variables:
    h3_10           (index) int64 91kB 622506101129773055 ... 622506153063612415
    height          (index) float64 91kB 221.2 237.3 242.0 ... 191.6 189.6 195.2
    sig0            (index) float64 91kB 25.33 20.98 36.63 ... 31.15 49.45 23.95
    classification  (index) float64 91kB 3.0 6.0 5.25 3.0 ... 3.0 3.0 3.0 3.0
    geoid           (index) float64 91kB 49.39 49.39 49.39 ... 49.33 49.33 49.33
    wse             (index) float64 91kB 171.8 187.9 192.6 ... 142.2 140.3 145.9
    geometry        (index) geometry 91kB <class 'xarray.core.extension_array...

xarray.Dataset

• Dimensions:
  • index: 11333
• Coordinates: (1)
  • index
    (index)
    int64
    0 1 2 3 ... 11329 11330 11331 11332

    array([    0,     1,     2, ..., 11330, 11331, 11332])

• Data variables: (7)
  • h3_10
    (index)
    int64
    622506101129773055 ... 622506153...

    array([622506101129773055, 622506101136523263, 622506101136654335, ...,
           622506153063350271, 622506153063383039, 622506153063612415])

  • height
    (index)
    float64
    221.2 237.3 242.0 ... 189.6 195.2

    array([221.17738342, 237.29124451, 241.98103714, ..., 191.55817159,
           189.62651571, 195.24961853])

  • sig0
    (index)
    float64
    25.33 20.98 36.63 ... 49.45 23.95

    array([25.32897758, 20.98498726, 36.63452053, ..., 31.14721298,
           49.45066833, 23.94847202])

  • classification
    (index)
    float64
    3.0 6.0 5.25 3.0 ... 3.0 3.0 3.0

    array([3.  , 6.  , 5.25, ..., 3.  , 3.  , 3.  ])

  • geoid
    (index)
    float64
    49.39 49.39 49.39 ... 49.33 49.33

    array([49.39186478, 49.39219666, 49.39214802, ..., 49.33113607,
           49.32905451, 49.33250427])

  • wse
    (index)
    float64
    171.8 187.9 192.6 ... 140.3 145.9

    array([171.78552246, 187.89904785, 192.58889008, ..., 142.22703552,
           140.29746501, 145.91711426])

  • geometry
    (index)
    geometry
    <class 'xarray.core.extension_ar...

    <class 'xarray.core.extension_array.PandasExtensionArray'>(array=<GeometryArray>
    [<POLYGON ((1.51 43.611, 1.51 43.611, 1.509 43.611, 1.51 43.61, 1.511 43.61, ...>,
     <POLYGON ((1.51 43.618, 1.51 43.618, 1.509 43.617, 1.51 43.616, 1.511 43.617...>,
     <POLYGON ((1.51 43.619, 1.509 43.619, 1.509 43.618, 1.51 43.618, 1.51 43.618...>,
     <POLYGON ((1.51 43.616, 1.509 43.616, 1.509 43.615, 1.51 43.614, 1.51 43.615...>,
     <POLYGON ((1.51 43.615, 1.51 43.614, 1.509 43.614, 1.51 43.613, 1.511 43.613...>,
     <POLYGON ((1.509 43.617, 1.509 43.617, 1.508 43.616, 1.509 43.616, 1.51 43.6...>,
     <POLYGON ((1.509 43.621, 1.508 43.621, 1.508 43.62, 1.509 43.62, 1.509 43.62...>,
     <POLYGON ((1.51 43.621, 1.51 43.621, 1.509 43.62, 1.51 43.62, 1.511 43.62, 1...>,
     <POLYGON ((1.508 43.62, 1.507 43.62, 1.507 43.619, 1.508 43.619, 1.508 43.61...>,
     <POLYGON ((1.51 43.622, 1.509 43.622, 1.509 43.621, 1.51 43.621, 1.51 43.621...>,
     ...
     <POLYGON ((1.374 43.663, 1.373 43.663, 1.373 43.663, 1.374 43.662, 1.374 43....>,
     <POLYGON ((1.371 43.663, 1.371 43.663, 1.37 43.662, 1.371 43.661, 1.372 43.6...>,
     <POLYGON ((1.377 43.663, 1.376 43.663, 1.376 43.662, 1.377 43.662, 1.377 43....>,
     <POLYGON ((1.376 43.664, 1.376 43.664, 1.375 43.663, 1.376 43.663, 1.377 43....>,
     <POLYGON ((1.378 43.664, 1.377 43.664, 1.377 43.663, 1.378 43.662, 1.378 43....>,
     <POLYGON ((1.375 43.663, 1.375 43.663, 1.374 43.662, 1.375 43.662, 1.376 43....>,
     <POLYGON ((1.375 43.667, 1.374 43.667, 1.374 43.667, 1.375 43.666, 1.375 43....>,
     <POLYGON ((1.376 43.667, 1.376 43.667, 1.375 43.666, 1.376 43.666, 1.377 43....>,
     <POLYGON ((1.375 43.669, 1.375 43.669, 1.374 43.669, 1.375 43.668, 1.376 43....>,
     <POLYGON ((1.371 43.664, 1.37 43.664, 1.37 43.663, 1.371 43.663, 1.371 43.66...>]
    Length: 11333, dtype: geometry)

• Indexes: (1)
  • index
    PandasIndex

    PandasIndex(RangeIndex(start=0, stop=11333, step=1, name='index'))

• Attributes: (0)

Display on folium map

Though not the most straightforward and efficient compared to QGIS or other python solutions, this would allow you to dispaly your data in leaflet starting from a GeoDataFrame

pixc_read.layers

['20230407_483_16_78L_10__h3', '20230406_482_16_78L_10__h3']

import folium
import branca.colormap as cmp

layer_h3 = pixc_read.read_single_layer('20230407_483_16_78L_10__h3')

# creating a colormap
linear = cmp.LinearColormap(
    ['blue', 'purple', 'orange', 'yellow'],
    vmin=layer_h3['wse'].min(),  # minimum wse value
    vmax=layer_h3['wse'].max(),  # maximum wse value
    caption='Water Surface Elevation (m)' #Caption for Color scale or Legend
)
# Initiating map
m = folium.Map([43.6, 1.43], zoom_start=12, tiles="cartodbpositron")

folium.GeoJson(
    layer_h3,
    style_function = lambda row:  {
        'fillColor': linear(row['properties']["wse"]),
        'weight': 0,          #how thick the border has to be
        'fillOpacity': 1
    },
).add_to(m)
linear.add_to(m)   #adds colorscale and legend
m

Make this Notebook Trusted to load map: File -> Trust Notebook

Enjoy !