Export to file#
A need for exporting dfsu files to other formats, such as raster formats (NetCDF, GeoTIFF) or vector formats (ESRI Shape, GeoJSON) are not uncommon.
There are usually very specific requirements about the exported data, naming, metadata, resolution and so on, which is why we don’t provide out of the box tools to convert files back and forth from a lot of formats, but instead provide the necessary building blocks and some examples to use as inspiration.
Examples#
Export to shapefile#
Exporting to GIS formats such as ESRI Shapefile or GeoTiff requires some libraries which are tricky to install on Windows.
One way is to use Conda-forge
Download the Miniconda installer
Install GeoPandas:
conda install --channel conda-forge geopandas
An alternative is to use to run python under Windows Subsystem for Linux
import mikeio
import geopandas as gpd
dfs = mikeio.open("data/oresundHD_run1.dfsu")
dfs
Dfsu2D
number of elements: 3612
number of nodes: 2046
projection: UTM-33
items:
0: Surface elevation <Surface Elevation> (meter)
1: Total water depth <Water Depth> (meter)
2: U velocity <u velocity component> (meter per sec)
3: V velocity <v velocity component> (meter per sec)
time: 5 steps with dt=86400.0s
2018-03-07 00:00:00 -- 2018-03-11 00:00:00
ds = dfs.read(items="Surface elevation")
ds
<mikeio.Dataset>
dims: (time:5, element:3612)
time: 2018-03-07 00:00:00 - 2018-03-11 00:00:00 (5 records)
geometry: Dfsu2D (3612 elements, 2046 nodes)
items:
0: Surface elevation <Surface Elevation> (meter)
maxwl = ds["Surface elevation"].max(axis=0).values
maxwl.shape
(3612,)
shp = dfs.geometry.to_shapely()
type(shp)
shapely.geometry.multipolygon.MultiPolygon
poly_list = [p for p in shp.geoms]
gdf = gpd.GeoDataFrame({'waterlevel': maxwl},geometry=poly_list, crs=dfs.projection_string)
---------------------------------------------------------------------------
CRSError Traceback (most recent call last)
Cell In[6], line 1
----> 1 gdf = gpd.GeoDataFrame({'waterlevel': maxwl},geometry=poly_list, crs=dfs.projection_string)
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/geopandas/geodataframe.py:192, in GeoDataFrame.__init__(self, data, geometry, crs, *args, **kwargs)
184 if (
185 hasattr(geometry, "crs")
186 and geometry.crs
187 and crs
188 and not geometry.crs == crs
189 ):
190 raise ValueError(crs_mismatch_error)
--> 192 self.set_geometry(geometry, inplace=True, crs=crs)
194 if geometry is None and crs:
195 raise ValueError(
196 "Assigning CRS to a GeoDataFrame without a geometry column is not "
197 "supported. Supply geometry using the 'geometry=' keyword argument, "
198 "or by providing a DataFrame with column name 'geometry'",
199 )
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/geopandas/geodataframe.py:349, in GeoDataFrame.set_geometry(self, col, drop, inplace, crs)
346 level.crs = crs
348 # Check that we are using a listlike of geometries
--> 349 level = _ensure_geometry(level, crs=crs)
350 # update _geometry_column_name prior to assignment
351 # to avoid default is None warning
352 frame._geometry_column_name = geo_column_name
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/geopandas/geodataframe.py:60, in _ensure_geometry(data, crs)
58 return GeoSeries(out, index=data.index, name=data.name)
59 else:
---> 60 out = from_shapely(data, crs=crs)
61 return out
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/geopandas/array.py:159, in from_shapely(data, crs)
143 def from_shapely(data, crs=None):
144 """
145 Convert a list or array of shapely objects to a GeometryArray.
146
(...)
157
158 """
--> 159 return GeometryArray(vectorized.from_shapely(data), crs=crs)
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/geopandas/array.py:294, in GeometryArray.__init__(self, data, crs)
291 self._data = data
293 self._crs = None
--> 294 self.crs = crs
295 self._sindex = None
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/geopandas/array.py:360, in GeometryArray.crs(self, value)
357 @crs.setter
358 def crs(self, value):
359 """Sets the value of the crs"""
--> 360 self._crs = None if not value else CRS.from_user_input(value)
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/pyproj/crs/crs.py:501, in CRS.from_user_input(cls, value, **kwargs)
499 if isinstance(value, cls):
500 return value
--> 501 return cls(value, **kwargs)
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/pyproj/crs/crs.py:348, in CRS.__init__(self, projparams, **kwargs)
346 self._local.crs = projparams
347 else:
--> 348 self._local.crs = _CRS(self.srs)
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/pyproj/_crs.pyx:2378, in pyproj._crs._CRS.__init__()
CRSError: Invalid projection: UTM-33: (Internal Proj Error: proj_create: several objects matching this name: Schwarzeck / UTM zone 33S, RGRDC 2005 / UTM zone 33S, Malongo 1987 / UTM zone 33S, ETRS89/DREF91/2016 / UTM zone 33N (N-zE), ETRS89/DREF91/2016 / UTM zone 33N (zE-N))
Ouch… The short and smart projection string “UTM-33” is apparently not understood by GeoPandas. Better look it up on the web https://spatialreference.org/ref/epsg/32433/
gdf = gpd.GeoDataFrame({'waterlevel': maxwl},geometry=poly_list, crs="EPSG:32433")
gdf.head()
| waterlevel | geometry | |
|---|---|---|
| 0 | 0.321391 | POLYGON ((353122.611 6199631.339, 354250.947 6... |
| 1 | 0.254802 | POLYGON ((356811.295 6165979.110, 356926.413 6... |
| 2 | 0.210542 | POLYGON ((334193.183 6143454.831, 332951.036 6... |
| 3 | 0.213418 | POLYGON ((330871.279 6142677.275, 331831.835 6... |
| 4 | 0.208291 | POLYGON ((334327.062 6140906.232, 335381.131 6... |
Export as ESRI Shapefile
gdf.to_file("waterlevel_utm.shp")
Which can be used together with other data sources and customized in QGIS
Or GeoJSON, which some might prefer…
gdf.to_file("waterlevel_utm.json")
