Extractor Tutorial

Real Data

The purpose of this tutorial is to show how to use the classes defined in xmsextractor to extract data at given locations on an XmUGrid with scalar values, using outputs from a numerical model.

First, we will load a UGrid containing elevation values. The data comes from a Sedimentation and River Hydraulics – Two-Dimensional (SRH-2D) model.

In [1]:
import numpy as np
import geoviews as gv
import holoviews as hv
import pandas as pd
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import panel as pn

from shapely.geometry import Polygon, shape
from holoviews.streams import Params

import param
import hvplot.pandas
import hvplot.xarray

import xmsextractor
import xmsgrid

import ugrid_tools
import extractor_tools
import shapefile

hv.extension('bokeh')

Display area of interest

In [2]:
srs = ccrs.LambertConformal(central_longitude=-98, central_latitude=35.0, false_easting=600000,
                           standard_parallels=(35.5667, 36.7667))
boundary_polygon = gv.Shape.from_shapefile("boundary_polys.shp", crs=srs).options(color="red", alpha=0.25)
basemap = gv.tile_sources.ESRI()

%opts Shape Polygons Points [height=700 width=700]
basemap * boundary_polygon
Out[2]:

Create a UGrid

In [3]:
extractor_example_data = extractor_tools.ExtractorExample()
ugrid = xmsgrid.ugrid.UGrid(extractor_example_data.points, extractor_example_data.cells)
View the UGrid
In [4]:
ugrid_cells = []

for i in range(0, ugrid.get_cell_count()):
    ugrid_cells.extend(ugrid.get_cell_cellstream(i)[1][2:])

polygons, lines, points = ugrid_tools.create_hv_objects_for_image(ugrid, _projection=srs)
basemap * polygons * lines * points
Out[4]:

Create an Extractor

In [5]:
extractor = xmsextractor.extractor.UGrid2dPolylineDataExtractor(ugrid, xmsextractor.extractor.data_location_enum.LOC_POINTS)

Setup Polyline

In [6]:
polyline = [(657675.3, 110106.7, 275.0), (657716.7, 110192.5, 275.0)]
extractor.set_polyline(polyline)
In [7]:
extracted_locations = extractor.get_extract_locations()

Set the point scalars for the water depth.

In [8]:
extractor.set_grid_scalars(extractor_example_data.wdepth_m, [], xmsextractor.extractor.data_location_enum.LOC_CELLS)
print(ugrid.get_locations())

[[6.58306098e+05 1.09672020e+05 2.75000000e+02]
 [6.58307941e+05 1.09679038e+05 2.75000000e+02]
 [6.58323570e+05 1.09697156e+05 2.75000000e+02]
 ...
 [6.58940545e+05 1.10397323e+05 2.70000000e+02]
 [6.58927408e+05 1.10385861e+05 2.70000000e+02]
 [6.58914221e+05 1.10373607e+05 2.70000000e+02]]

Extract the data.

In [9]:
extracted_data = extractor.extract_data()
print(str(["{0:0.4f}".format(i) for i in extracted_data]).replace("'", ""))

[0.4364, 0.4364, 0.7657, 0.9219, 1.3677, 1.5640, 1.7878, 1.8955, 1.9747, 2.0098, 2.0801, 2.1397, 2.1662, 2.2067, 2.2315, 2.2453, 2.2668, 2.2565, 2.2500, 2.2132, 2.2099, 2.1895, 1.8264, 1.4297, 1.4297]

View the Data

In [10]:
extract_scatter_vals = np.array([
    (extracted_locations[i][0], extracted_locations[i][1], extracted_data[i])
    for i in range(0, len(extracted_data))
])

min_x = min([x[0] for x in polyline])
min_y = min([y[1] for y in polyline])
max_x = max([x[0] for x in polyline])
max_y = max([y[1] for y in polyline])

padding = (max((max_x-min_x), (max_y-min_y))) * 0.15

#hv_extract_values = gv.Scatter(extract_scatter_vals, crs=srs, vdims=['y','z'])
hv_polyline = gv.Path([polyline,], crs=srs).options(line_color='red')
%opts Polygons Points Path [height=500 width=500]
%opts Scatter [color_index=2 colorbar=True tools=['hover']] (size=10)
#polygons, lines, points = ugrid_tools.create_hv_objects_for_image(ugrid, _projection=srs)
polygons = polygons.redim.range(
    x=(min_x-padding, max_x+padding),
    y=(min_y-padding, max_y+padding)
)
In [11]:
# draw the geoviews map
basemap * polygons * lines * points * hv_polyline * hv_extract_values

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-11-e3c647bcdc8c> in <module>
      1 # draw the geoviews map
----> 2 basemap * polygons * lines * points * hv_polyline * hv_extract_values

NameError: name 'hv_extract_values' is not defined
In [12]:
pxn = extracted_locations[-1][0]
px0 = extracted_locations[0][0]
distance = pxn - px0
hv_point_locations = [(extracted_locations[i][0] - px0, extracted_data[i]) for i in range(0, len(extracted_data))]
height = max([y[1] for y in hv_point_locations])

x_range = (min([x[0] for x in hv_point_locations]) - (distance*0.1), max([x[0] for x in hv_point_locations]) + (distance*0.1))
y_range = (min([y[1] for y in hv_point_locations]) - (height*0.1),   max([y[1] for y in hv_point_locations]) + (height*0.1))

%opts Points [colorbar=True tools=['hover']] (size=4 color='red')
graph_points = hv.Points(hv_point_locations)
graph_line = hv.Path([hv_point_locations,],)
graph_line * graph_points.redim.range(x=x_range, y=y_range)
Out[12]: