xmsextractor

Introduction

The purpose of this tutorial is to provide explanation on how to use the classes defined in xmsextractor to extract data at given locations on an unstructured grid (xmsgrid) with scalar values.

The motivation for this library is a tool to interrogate results from numerical models. Numerical models use either structured or unstructured grids to perform calculations. The locations where the numerical model performs calculations do not usually coincide with points of interest. This library extracts data values from arbitrary locations in an unstructured grid (or structured grid); also, the library includes functionality to extract values along a line.

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
from holoviews.streams import Params

import xmsextractor
import xmsgrid

import ugrid_tools

hv.extension('bokeh')

Example - Simple Location Extractor

This first example shows how to extract data from a simple XmUGrid. A picture of the example is shown below. Notice that the XmUGrid is a simple square from (0,0) to (1,1) split into two triangle cells. The scalar data is mapped to the cells (this usually means that the numerical model has computed a value at the cell center).

If the scalar data were mapped to the points then that would mean that the numerical model computed values at the point locations.

In [2]:
# create a simple xmsgrid
pts = [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)]
cells = [
    xmsgrid.ugrid.UGrid.TRIANGLE, 3, 0, 1, 2,
    xmsgrid.ugrid.UGrid.TRIANGLE, 3, 2, 3, 0,
]
ugrid = xmsgrid.ugrid.UGrid(pts, cells)
In [3]:
# draw the grid
polygons, lines, points = ugrid_tools.create_hv_objects_for_image(ugrid) 
%opts Polygons Points [width=700, height=700]
polygons.options(color="red", alpha=0.25) * lines * points.options(size=3)
Out[3]:

Data extraction at locations

The basic steps to extract interpolated locations from an XmUGrid and scalar values is as follows:

  1. Create an extractor for an existing XmUGrid.
In [4]:
extractor = xmsextractor.extractor.UGrid2dDataExtractor(ugrid)
  1. Set scalar and activity values.

Usually these values would come from a numerical model. Activity values refers to output from a numerical model where certain areas in the grid do not have computed values. For a 2D simulation of the flow of water the "inactive" areas could be areas that are dry.

In [5]:
cellScalars = (1, 2)
extractor.set_grid_cell_scalars(cellScalars, (), xmsextractor.extractor.data_location_enum.LOC_CELLS)
  1. Set extract locations.

These are the locations where we want to know the value computed by the numerical model. The value will be interpolated based on the surrounding values computed by the numerical model.

In [6]:
extractor_locations = [
    (0.0, 0.0, 0.0),
    (0.25, 0.75, 100.0),
    (0.5, 0.5, 0.0),
    (0.75, 0.25, -100.0),
    (-0.1, -0.1, 0.0),
    (0, 1, 0),
    (1, 1, 0),
    (1, 0, 0),
    (0.25, 0.4, 0),
    (0.75, 0.6, 0),
]

# draw the locations
extractor_location_points = hv.Points(extractor_locations).options(color='red', size=10)
%opts Polygons Points [width=700, height=700]
polygons * lines * points * extractor_location_points
Out[6]:
In [7]:
extractor.set_extract_locations(extractor_locations)
  1. Extract the data. Values that are outside of the UGrid are returned as nan by default.
In [9]:
extraction_values = extractor.extract_data()
for i in range(0, len(extractor_locations)):
    print(f"Extraction value at {extractor_locations[i]}: {extraction_values[i]}")

Extraction value at (0.0, 0.0, 0.0): 1.5
Extraction value at (0.25, 0.75, 100.0): 2.0
Extraction value at (0.5, 0.5, 0.0): 1.5
Extraction value at (0.75, 0.25, -100.0): 1.0
Extraction value at (-0.1, -0.1, 0.0): nan
Extraction value at (0, 1, 0): 2.0
Extraction value at (1, 1, 0): 1.5
Extraction value at (1, 0, 0): 1.0
Extraction value at (0.25, 0.4, 0): 1.725000023841858
Extraction value at (0.75, 0.6, 0): 1.274999976158142
View the Data
In [10]:
extract_scatter_vals = np.array([
    (extractor_locations[i][0], 
     extractor_locations[i][1], 
     extraction_values[i])
    for i in range(0, len(extraction_values))
])
hv_extract_values = hv.Scatter(extract_scatter_vals, vdims=['y', 'z'])
polygons, lines, points = ugrid_tools.create_hv_objects_for_image(ugrid) 
%opts Polygons Points [width=700, height=700 tools=['hover']]
%opts Scatter [color_index=2 colorbar=True tools=['hover']] (size=10 cmap="cool")
polygons * lines * points.options(size=3) * hv_extract_values
Out[10]:

Example - Transient Data Location Extractor

This example shows how to use a location extractor on transient data.

This example shows how to extract data from a XmUGrid that has transient scalar data. A picture of the example is shown below. Notice that the UGrid is a 2x3 structured grid with quadrillateral cells. There are two time steps with scalar data mapped to the points.

1. Create a UGrid for the Extractor

In [11]:
transient_example_points = [(288050, 3907770, 0), (294050, 3907770, 0), (300050, 3907770, 0),
          (306050, 3907770, 0), (288050, 3901770, 0), (294050, 3901770, 0),
          (300050, 3901770, 0), (306050, 3901770, 0), (288050, 3895770, 0),
          (294050, 3895770, 0), (300050, 3895770, 0), (306050, 3895770, 0)
]
transient_example_cells = [xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 0, 4, 5, 1,
         xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 1, 5, 6, 2,
         xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 2, 6, 7, 3,
         xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 4, 8, 9, 5,
         xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 5, 9, 10, 6,
         xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 6, 10, 11, 7
]
transient_example_ugrid = xmsgrid.ugrid.UGrid(transient_example_points, transient_example_cells)

polygons, lines, points = ugrid_tools.create_hv_objects_for_image(transient_example_ugrid) 
%opts Polygons Points [width=700, height=700]
polygons * lines * points.options(size=3)
Out[11]:

2. Create the Extractor

In [12]:
transient_example_extractor = xmsextractor.extractor.UGrid2dDataExtractor(transient_example_ugrid)

3. Set extract locations.

In [13]:
transient_example_extract_locations = [(289780, 3906220, 0),
                     (293780, 3899460, 0),
                     (298900, 3900780, 0),
                     (301170, 3904960, 0),
                     (296330, 3906180, 0),
                     (307395, 3901463, 0)
]
transient_example_extractor.set_extract_locations(transient_example_extract_locations)
Display the Extract Locations
In [14]:
transient_example_extractor_location_points = hv.Points(transient_example_extract_locations).options(color='red', size=10)
%opts Polygons Points [width=700, height=700]
polygons * lines * points * transient_example_extractor_location_points
Out[14]:

4. Set the point scalars for the first time step.

In [16]:
transient_example_point_scalars = [730.787, 1214.54, 1057.145, 629.2069, 351.1153, 631.6649, 1244.366,
                                   449.9133, 64.04247, 240.9716, 680.0491, 294.9547]
transient_example_extractor.set_grid_point_scalars(
    transient_example_point_scalars, [], xmsextractor.extractor.data_location_enum.LOC_CELLS
)

5. Extract the data.

In [17]:
transient_example_extraction_values = transient_example_extractor.extract_data()
for i in range(0, len(transient_example_extract_locations)):
    print(f"Extraction value at {transient_example_extract_locations[i]}: {transient_example_extraction_values[i]}")

Extraction value at (289780, 3906220, 0): 719.6930541992188
Extraction value at (293780, 3899460, 0): 468.6232604980469
Extraction value at (298900, 3900780, 0): 1033.8193359375
Extraction value at (301170, 3904960, 0): 996.5289916992188
Extraction value at (296330, 3906180, 0): 1204.343505859375
Extraction value at (307395, 3901463, 0): nan

6. Continue using steps 4 and 5 for remaining time steps.

In [18]:
transient_example_point_scalars = [float('nan'), 1220.5, 1057.1, 613.2, 380.1, 625.6, 722.2, 449.9, 51.0, 240.9, 609.0, 294.9]
transient_example_cell_activity = [True] * transient_example_ugrid.get_cell_count()
transient_example_cell_activity[0] = False
transient_example_extractor.set_grid_point_scalars(
    transient_example_point_scalars, 
    transient_example_cell_activity, 
    xmsextractor.extractor.data_location_enum.LOC_CELLS
)

transient_example_extraction_values = transient_example_extractor.extract_data()
for i in range(0, len(transient_example_extract_locations)):
    print(f"Extraction value at {transient_example_extract_locations[i]}: {transient_example_extraction_values[i]}")

Extraction value at (289780, 3906220, 0): nan
Extraction value at (293780, 3899460, 0): 466.4429931640625
Extraction value at (298900, 3900780, 0): 685.0070190429688
Extraction value at (301170, 3904960, 0): 849.4258422851562
Extraction value at (296330, 3906180, 0): 1069.6595458984375
Extraction value at (307395, 3901463, 0): nan

View the Data

In [19]:
transient_extract_scatter_vals = np.array([
    (transient_example_extract_locations[i][0], 
     transient_example_extract_locations[i][1], 
     transient_example_extraction_values[i])
    for i in range(0, len(transient_example_extraction_values))
])
hv_transient_extract_values = hv.Scatter(transient_extract_scatter_vals, vdims=['y', 'z'])
polygons, lines, points = ugrid_tools.create_hv_objects_for_image(transient_example_ugrid) 
%opts Polygons Points [width=700, height=700 tools=['hover']]
%opts Scatter [color_index=2 colorbar=True tools=['hover']] (size=10)
polygons * lines * points.options(size=3) * hv_transient_extract_values
Out[19]:

Example - Transient Polyline Extractor

This example shows how to use a polyline extractor on transient data.

This example shows how to extract data along a polyline from a XmUGrid that has transient scalar data. A picture of the example is shown below. Notice that the UGrid is a 2x3 structured grid with quadrillateral cells. There are two time steps with scalar data mapped to the points.

The steps to extract interpolated values along a polyline for transient scalar values include:

1. Build the UGrid and Extractor Objects

In [20]:
polyline_example_points = [(288050, 3907770, 0), (294050, 3907770, 0), (300050, 3907770, 0),
          (306050, 3907770, 0), (288050, 3901770, 0), (294050, 3901770, 0),
          (300050, 3901770, 0), (306050, 3901770, 0), (288050, 3895770, 0),
          (294050, 3895770, 0), (300050, 3895770, 0), (306050, 3895770, 0)
]
polyline_example_cells = [xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 0, 4, 5, 1,
         xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 1, 5, 6, 2,
         xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 2, 6, 7, 3,
         xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 4, 8, 9, 5,
         xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 5, 9, 10, 6,
         xmsgrid.ugrid.UGrid.ugrid_celltype_enum.QUAD, 4, 6, 10, 11, 7
]
polyline_example_ugrid = xmsgrid.ugrid.UGrid(polyline_example_points, polyline_example_cells)

# Step 1. Create an extractor for an XmUGrid giving the mapped location of the scalar values
polyline_example_extractor = xmsextractor.extractor.UGrid2dPolylineDataExtractor(
    polyline_example_ugrid, xmsextractor.extractor.data_location_enum.LOC_POINTS
)

polygons, lines, points = ugrid_tools.create_hv_objects_for_image(polyline_example_ugrid) 
%opts Polygons Points [width=700, height=700]
polygons * lines * points.options(size=3)
Out[20]:

2. Set the polyline

In [22]:
# Step 3. Set the polyline to be extracted along.
polyline = [(290764, 3895106, 0), (291122, 3909108, 0),
            (302772, 3909130, 0), (302794, 3895775, 0)
]
polyline_example_extractor.set_polyline(polyline)

hv_polyline = hv.Path([polyline,]).options(line_color='red')
hv_polyline_points = hv.Points(polyline).options(size=6)
polygons, lines, points = ugrid_tools.create_hv_objects_for_image(polyline_example_ugrid) 
%opts Polygons Points [width=700, height=700]
polygons * lines * points.options(size=3) * hv_polyline * hv_polyline_points
Out[22]:
Display extraction locations generated by the polyline
In [23]:
polyline_example_extracted_locations = polyline_example_extractor.get_extract_locations()

hv_polyline = hv.Path([polyline,]).options(line_color='red')
hv_polyline_points = hv.Points(polyline_example_extracted_locations).options(size=10)
polygons, lines, points = ugrid_tools.create_hv_objects_for_image(polyline_example_ugrid) 
%opts Polygons Points [width=700, height=700]
polygons * lines * points.options(size=3) * hv_polyline * hv_polyline_points
Out[23]:

3. Set the point scalars for the first time step.

In [24]:
point_scalars = [730.787, 1214.54, 1057.145, 629.2069, 351.1153, 631.6649,
                 1244.366, 449.9133, 64.04247, 240.9716, 680.0491, 294.9547
]
polyline_example_extractor.set_grid_scalars(point_scalars, [], xmsextractor.extractor.data_location_enum.LOC_CELLS)

4. Extract the data

In [25]:
polyline_example_extraction_values = polyline_example_extractor.extract_data()
for i in range(0, len(polyline_example_extracted_locations)):
    print(f"Extraction value at {polyline_example_extracted_locations[i]}: {polyline_example_extraction_values[i]}")

Extraction value at [ 290764. 3895106.       0.]: nan
Extraction value at [ 290780.97700329 3895770.               0.        ]: 144.57403564453125
Extraction value at [ 290862.47493036 3898957.52506964       0.        ]: 299.48590087890625
Extraction value at [ 290934.38365948 3901770.               0.        ]: 485.9841003417969
Extraction value at [ 291012.05710306 3904807.94289694       0.        ]: 681.852783203125
Extraction value at [ 291087.79031567 3907770.               0.        ]: 975.7103881835938
Extraction value at [ 291122. 3909108.       0.]: nan
Extraction value at [ 302772. 3909130.       0.]: nan
Extraction value at [ 302774.24035942 3907770.               0.        ]: 862.843994140625
Extraction value at [ 302778.73546839 3905041.26453161       0.        ]: 780.9825439453125
Extraction value at [ 302784.12429802 3901770.               0.        ]: 882.3439331054688
Extraction value at [ 302788.63571589 3899031.36428411       0.        ]: 811.0173950195312
Extraction value at [ 302794. 3895775.       0.]: 504.4028625488281

5. Continue using steps 3 and 4 for remaining time steps.

In [26]:
point_scalars = [float('nan'), 1220.5, 1057.1, 613.2, 380.1, 625.6, 722.2, 449.9, 51.0, 240.9,
                 609.0, 294.9]
polyline_example_extractor.set_grid_scalars(point_scalars, [], xmsextractor.extractor.data_location_enum.LOC_CELLS)
extract_vals = polyline_example_extractor.extract_data()

View the Data

In [27]:
extract_scatter_vals = np.array([
    (polyline_example_extracted_locations[i][0], polyline_example_extracted_locations[i][1], extract_vals[i])
    for i in range(0, len(extract_vals))
])
hv_extract_values = hv.Scatter(extract_scatter_vals, vdims=['y', 'z'])
polygons, lines, points = ugrid_tools.create_hv_objects_for_image(polyline_example_ugrid) 
%opts Polygons Points [width=700, height=700 tools=['hover']]
%opts Scatter [color_index=2 colorbar=True tools=['hover']] (size=10)
polygons * lines * points.options(size=3) * hv_polyline * hv_extract_values
Out[27]: