Feature/ugrid auto mesh

Dockerfile

@@ -5,9 +5,10 @@ ENV LANG=C.UTF-8 LC_ALL=C.UTF-8
 RUN \
     apt-get update --fix-missing && \
     apt-get install -y ffmpeg wget unzip libglu1-mesa-dev gcc nfs-common && \
-	&& rm -rf /var/lib/apt/lists/*
+    rm -rf /var/lib/apt/lists/*
 
-RUN conda install -c conda-forge libgdal gdal
+RUN conda config --add channels conda-forge
+RUN conda install libgdal gdal vtk mayavi netcdf4 numpy scipy scikit-image matplotlib pandas rasterio Shapely pillow cython
 
 # install flowmap
 COPY ./ app/

flowmap/cli.py

@@ -254,7 +254,6 @@ def export(dataset, dem, format, **kwargs):
         )
 
 
-
     klass = flowmap.formats.get_format(dataset, **kwargs)
     ds = klass(dataset, dem=dem, **kwargs)
     logger.info("exporting grid")

flowmap/dem.py

@@ -15,9 +15,10 @@ def read_dem(dem_filename):
         dem['band'] = band
 
     dem['transform'] = src.get_transform()
-    dem['affine'] = src.affine
+    # is the affine transformation
+    dem['affine'] = src.transform
     # pixel sizes
-    affine = src.affine
+    affine = dem['affine']
     dem['dxp'] = affine.a
     dem['dyp'] = -affine.e
     dem['width'] = src.width

flowmap/formats/ugrid.py

@@ -12,7 +12,9 @@
 
 import netCDF4
 import numpy as np
-import pyugrid
+# rename because I already named my variable ugrid
+from gridded.pyugrid import ugrid as pyugrid
+from gridded.pyugrid import read_netcdf
 import geojson
 
 # used for transforming into a vtk grid and for particles
@@ -23,6 +25,10 @@
 import rasterio.crs
 import shapely.geometry
 
+from shapely import speedups
+# TODO: check if needed in container
+speedups.disable()
+
 from .netcdf import NetCDF
 from .. import particles
 from .. import subgrid
@@ -53,27 +59,32 @@ class UGrid(NetCDF):
 
     def validate(self):
         """validate a file"""
-        valid = True
+        valid = False
 
         with netCDF4.Dataset(self.path) as ds:
-            variables = ds.variables.keys()
-        # for now we hardcode the filenames. This can be replaced by the pyugrid, once released
-        for var in ("mesh2d_face_nodes", "mesh2d_node_x", "mesh2d_node_x"):
-            if var not in variables:
-                logger.warn(
-                    "%s not found in variables of file %s",
-                    var,
-                    self.path
-                )
-                valid = False
-                return valid
+            conventions = getattr(ds, 'Conventions', '')
+            if 'ugrid' in conventions.lower():
+                valid = True
         return valid
 
+
+    @property
+    def mesh2d(self):
+        with netCDF4.Dataset(self.path) as ds:
+            mesh_names = read_netcdf.find_mesh_names(ds)
+            if not mesh_names:
+                raise ValueError('No meshes found in {}'.format(self.path))
+            for mesh_name in mesh_names:
+                var = ds.variables[mesh_name]
+                if var.topology_dimension == 2:
+                    return mesh_name
+        raise ValueError('No 2d mesh found in {}'.format(mesh_names))
+
     @property
     def ugrid(self):
         """Generate a ugrid grid from the input"""
         # TODO, lookup mesh name
-        ugrid = pyugrid.UGrid.from_ncfile(self.path, 'mesh2d')
+        ugrid = pyugrid.UGrid.from_ncfile(self.path, self.mesh2d)
 
         faces = np.ma.asanyarray(ugrid.faces)
 
@@ -105,16 +116,21 @@ def to_features(self):
         ugrid = self.ugrid
         faces = ugrid['faces']
 
+        # removed in geojson >= 2.5
         crs = geojson.crs.Named(
             properties={
                 "name": "urn:ogc:def:crs:EPSG::{:d}".format(self.src_epsg)
             }
         )
-        counts = (~faces.mask).sum(axis=1)
+        mask = np.ma.getmaskarray(faces)
+        counts = (~mask).sum(axis=1)
         face_coordinates = ugrid['face_coordinates']
         features= []
         for i, (face, count) in tqdm.tqdm(enumerate(zip(face_coordinates, counts)), desc='grid->features'):
-            poly = shapely.geometry.Polygon(face[:count].tolist())
+            coords = face[:count].tolist()
+            # close hull
+            coords.append(coords[0])
+            poly = shapely.geometry.Polygon(coords)
             geometry = poly.__geo_interface__
             geometry['crs'] = dict(crs)
             feature = geojson.Feature(id=i, geometry=geometry)
@@ -131,26 +147,43 @@ def to_polys(self):
                 "name": "urn:ogc:def:crs:EPSG::{:d}".format(self.src_epsg)
             }
         )
-        counts = (~faces.mask).sum(axis=1)
+
+        mask = np.ma.getmaskarray(faces)
+        counts = (~mask).sum(axis=1)
         face_coordinates = ugrid['face_coordinates']
         polys = []
         for i, (face, count) in tqdm.tqdm(enumerate(zip(face_coordinates, counts)), desc='grid->polys'):
-            poly = shapely.geometry.Polygon(face[:count].tolist())
+            coords = face[:count].tolist()
+            # close hull
+            coords.append(coords[0])
+            poly = shapely.geometry.Polygon(coords)
             polys.append(poly)
         return polys
 
-    def to_polydata(self):
+
+    def to_polydata(self, transform=False):
         """convert grid to a vtk polydata object"""
         ugrid = self.ugrid
 
         faces = ugrid['faces']
+
         points = ugrid['points']
+        # transform points
+        if transform:
+            points = np.array(
+                self.srs['src2utm'].TransformPoints(
+                    np.c_[points[:, 1], points[:, 0]]
+                )
+            )
+
+
 
         n_cells = faces.shape[0]
 
         cell_array = tvtk.CellArray()
 
-        counts = (~faces.mask).sum(axis=1)
+        mask = np.ma.getmaskarray(faces)
+        counts = (~mask).sum(axis=1)
         assert faces.min() >= 0, 'expected 0 based faces'
         cell_idx = np.c_[counts, faces.filled(-999)].ravel()
         cell_idx = cell_idx[cell_idx != -999]
@@ -177,9 +210,9 @@ def waterlevel(self, t):
         """lookup the waterlevel, depth and volume on timestep t"""
         # TODO: inspect mesh variable
         with netCDF4.Dataset(self.path) as ds:
-            s1 = ds.variables['mesh2d_s1'][t]
-            waterdepth = ds.variables['mesh2d_waterdepth'][t]
-            vol1 = ds.variables['mesh2d_vol1'][t]
+            s1 = ds.variables[self.mesh2d + '_s1'][t]
+            waterdepth = ds.variables[self.mesh2d + '_waterdepth'][t]
+            vol1 = ds.variables[self.mesh2d + '_vol1'][t]
         return dict(
             s1=s1,
             vol1=vol1,
@@ -191,8 +224,9 @@ def velocities(self, t):
         # TODO: inspect mesh variables
         with netCDF4.Dataset(self.path) as ds:
             # cumulative velocities
-            ucx = ds.variables['mesh2d_ucx'][t]
-            ucy = ds.variables['mesh2d_ucy'][t]
+            ucx = ds.variables[self.mesh2d + '_ucx'][t]
+            ucy = ds.variables[self.mesh2d + '_ucy'][t]
+
         return dict(
             ucx=ucx,
             ucy=ucy

requirements.txt

@@ -1,18 +1,7 @@
 Click==7.0
-mayavi==4.6.2
-netcdf4==1.4.2
-scikit-image==0.14.1
-numpy==1.15.4
-scipy==1.1.0
 tqdm==4.28.1
+# don't upgrade
 geojson==2.4.1
 mako==1.0.7
-matplotlib==3.0.2
-pandas==0.23.4
-rasterio==1.0.12
-geojson==2.4.1
-pyugrid==0.3.1
-Shapely==1.6.4
+https://github.com/NOAA-ORR-ERD/gridded/archive/v0.2.5.tar.gz
 simplejson==3.16.0
-cython==0.29.1
-pillow==5.3.0

requirements_dev.txt

@@ -1,10 +1,9 @@
-pip==18.1
+pip==20.2.3
 bumpversion==0.5.3
-wheel==0.30.0
+wheel==0.35.1
 watchdog==0.8.3
 flake8==3.5.0
 coverage==4.4.2
 Sphinx==1.6.5
 cryptography>=2
-PyYAML==3.13
-awscli==1.16.87
+awscli==1.18.154