Support for ADCIRC file [UGRID]

[marceloandrioni]

Hello, I am trying to load ADCIRC hydrodynamic model results in QGIS using "Add Mesh Layer" but it's not working.

After some small tweaks to make the files more UGRID (unstructured grid) compatible I was able to load them using EDAL and ncWMS as show in the image below, so I am quite sure the files are OK.

The files used in the test are here: ugrids.zip

Thank you.

[PeterPetrik]

@marceloandrioni hi. it looks like the elements variable is missing _FillValue attribute. I have a question: attached example is directly output of the solver, or is it modified version? If it is modified, can you please attach the original example?

[marceloandrioni]

Hello @PeterPetrik, I tried adding a _FillValue to the element variable but the same error (Invalid Data Source) happened.
Answering your question, in the first message I attached a slightly modified version of the raw output. I had to do theses changes so that the file could be read by EDAL/ncWMS. Basically, I changed the time:units and some attributes in the mesh_topology to follow CF/UGRID.

Here is an example of the direct output of ADCIRC. Note: this file is actually converted by a ADCIRC auxiliary script (https://adcirc.org/home/related-software/adcirc-utility-programs/) from the raw ascii output.

Here some other example hosted on a THREDDS by NOAA (but I think this one had some changes from the original format also).

[PeterPetrik]

could you please cut some timesteps from the direct output so it has < 4 MB ? I will have the fix ready today

[PeterPetrik]

obviously the original file is not-ugrid/cf conforming:
for beginning it misses: cf_role = mesh_topology attribute, etc. to define the mesh topology
also it looks like the variable/array names are bit non-standard, so I am wondering if it make sense to write custom parser for such non-conforming files.

[PeterPetrik]

@marceloandrioni if you know how to build you qgis with MDAL (you need to specify cmake's flag WITH_INTERNAL_MDAL=off and set your MDAL_PREFIX for custom mdal build), the fix from #156 should allow you to open your modified files

[marceloandrioni]

Hello @PeterPetrik, thank you for the lightning-fast response. Sorry for the delay on my part, but I was at a meeting without my computer.

Here are two small files with just two timesteps. One is the original netcdf output from ADCIRC (adcirc_original_2times.nc) and the other the modified version (adcirc_modified_2times.nc). I completely agree with you that writing a parser to deal with the original output of ADCIRC that doesn't follow the UGRID/CF conventions doesn't make sense.
But I think that being able to read the slightly modified version (that follows CF/UGRID) I created will be a really good improvement for the ADCIRC modeling community.

Just to leave a register for someone else dealing with this problem, here is the python script I made that puts ADCIRC results in CF/UGRID format.

To answer your last question, unfortnely I never built qgis from source, I just use the last version from the repos. So I think I will have to wait for the next release to acctualy use the improved version.

Thank you very much.

[BrianOBlanton]

Hi @PeterPetrik, I'm an ADCIRC user and in the dev group for ADCIRC. A colleague pointed me to MDAL/QGIS, which is very cool. I know very little about GIS but was able to easily get an ADCIRC output netCDF file into QGIS (3.12.2). This is a great resource for connecting ADCIRC (and other grid models) to the GIS community.

I'm curious about the issues people have noted re: CF-UGRID compliance. Older ADCIRC files may not be, but more recent files definitely should be. ADCIRC optionally outputs variables in netCDF format with CF-UGRID metadata.

I'd be very happy to help sort out any compliance issues if that would be useful.

Cheers,
Brian

[marceloandrioni]

Hi @BrianOBlanton, I was the one that reported the problems loading ADCIRC files into QGIS/MDAL. The ADCIRC runs were not done by me, so I am not sure about the version used. When you say recent files should be ok, do you know in which version of ADCIRC the changes in the NetCDF/UGRID attributes were made?

The output files I have (from a daily operational run) are the fort.63.nc (zeta) and fort.64.nc (uv), with the following attributes:

netcdf ADCIRC_BG_zeta_20200312 {
dimensions:
	time = UNLIMITED ; // (191 currently)
	node = 12769 ;
	nele = 23860 ;
	nvertex = 3 ;
	single = 1 ;
	nope = 1 ;
	max_nvdll = 31 ;
	neta = 31 ;
	nbou = 54 ;
	nvel = 1806 ;
	max_nvell = 932 ;
variables:
	double time(time) ;
		time:long_name = "model time" ;
		time:standard_name = "time" ;
		time:units = "12769" ;
	double x(node) ;
		x:long_name = "longitude" ;
		x:standard_name = "longitude" ;
		x:units = "degrees_east" ;
		x:positive = "east" ;
	double y(node) ;
		y:long_name = "latitude" ;
		y:standard_name = "latitude" ;
		y:units = "degrees_north" ;
		y:positive = "north" ;
	int element(nele, nvertex) ;
		element:long_name = "element" ;
		element:standard_name = "face_node_connectivity" ;
		element:units = "nondimensional" ;
		element:start_index = 1 ;
	double nvdll(nope) ;
		nvdll:long_name = "total number of nodes in each elevation specified & boundary segment" ;
		nvdll:units = "nondimensional" ;
	double nbdv(max_nvdll, nope) ;
		nbdv:long_name = "node numbers on each elevation specified boundary & segment" ;
		nbdv:units = "nondimensional" ;
	double nvell(nbou) ;
		nvell:long_name = "number of nodes in each normal flow specified boundary segment" ;
		nvell:units = "nondimensional" ;
	double ibtype(nbou) ;
		ibtype:long_name = "type of normal flow (discharge) boundary" ;
		ibtype:units = "nondimensional" ;
	double nbvv(max_nvell, nbou) ;
		nbvv:long_name = "node numbers on normal flow boundary segment" ;
		nbvv:units = "nondimensional" ;
	double depth(node) ;
		depth:long_name = "distance from geoid" ;
		depth:standard_name = "depth_below_geoid" ;
		depth:coordinates = "time y x" ;
		depth:location = "node" ;
		depth:mesh = "adcirc_mesh" ;
		depth:units = "m" ;
		depth:positive = "down" ;
	int adcirc_mesh(single) ;
		adcirc_mesh:long_name = "mesh topology" ;
		adcirc_mesh:standard_name = "mesh_topology" ;
		adcirc_mesh:dimension = 2 ;
		adcirc_mesh:node_coordinates = "x y" ;
		adcirc_mesh:face_node_connectivity = "element" ;
	double zeta(time, node) ;
		zeta:_FillValue = -99999. ;
		zeta:long_name = "water surface elevation above geoid" ;
		zeta:standard_name = "sea_surface_height_above_geoid" ;
		zeta:coordinates = "time y x" ;
		zeta:location = "node" ;
		zeta:mesh = "adcirc_mesh" ;
		zeta:units = "m" ;

// global attributes:
		:primitive_weighting_in_continuity_equation = "unitless" ;
		:sea_surface_height_above_geoid = "meters" ;
		:average_horizontal_eddy_viscosity_in_sea_water_wrt_depth = "(length**2)" ;
}

and

netcdf ADCIRC_BG_uv_20200312 {
dimensions:
	time = UNLIMITED ; // (191 currently)
	node = 12769 ;
	nele = 23860 ;
	nvertex = 3 ;
	single = 1 ;
	nope = 1 ;
	max_nvdll = 31 ;
	neta = 31 ;
	nbou = 54 ;
	nvel = 1806 ;
	max_nvell = 932 ;
variables:
	double time(time) ;
		time:long_name = "model time" ;
		time:standard_name = "time" ;
		time:units = "12769" ;
	double x(node) ;
		x:long_name = "longitude" ;
		x:standard_name = "longitude" ;
		x:units = "degrees_east" ;
		x:positive = "east" ;
	double y(node) ;
		y:long_name = "latitude" ;
		y:standard_name = "latitude" ;
		y:units = "degrees_north" ;
		y:positive = "north" ;
	int element(nele, nvertex) ;
		element:long_name = "element" ;
		element:standard_name = "face_node_connectivity" ;
		element:units = "nondimensional" ;
		element:start_index = 1 ;
	double nvdll(nope) ;
		nvdll:long_name = "total number of nodes in each elevation specified & boundary segment" ;
		nvdll:units = "nondimensional" ;
	double nbdv(max_nvdll, nope) ;
		nbdv:long_name = "node numbers on each elevation specified boundary & segment" ;
		nbdv:units = "nondimensional" ;
	double nvell(nbou) ;
		nvell:long_name = "number of nodes in each normal flow specified boundary segment" ;
		nvell:units = "nondimensional" ;
	double ibtype(nbou) ;
		ibtype:long_name = "type of normal flow (discharge) boundary" ;
		ibtype:units = "nondimensional" ;
	double nbvv(max_nvell, nbou) ;
		nbvv:long_name = "node numbers on normal flow boundary segment" ;
		nbvv:units = "nondimensional" ;
	double depth(node) ;
		depth:long_name = "distance from geoid" ;
		depth:standard_name = "depth_below_geoid" ;
		depth:coordinates = "time y x" ;
		depth:location = "node" ;
		depth:mesh = "adcirc_mesh" ;
		depth:units = "m" ;
	int adcirc_mesh(single) ;
		adcirc_mesh:long_name = "mesh topology" ;
		adcirc_mesh:standard_name = "mesh_topology" ;
		adcirc_mesh:dimension = 2 ;
		adcirc_mesh:node_coordinates = "x y" ;
		adcirc_mesh:face_node_connectivity = "element" ;
	double u-vel(time, node) ;
		u-vel:_FillValue = -99999. ;
		u-vel:long_name = "water column vertically averaged east/west velocity" ;
		u-vel:standard_name = "barotropic_eastward_sea_water_velocity" ;
		u-vel:coordinates = "time y x" ;
		u-vel:location = "node" ;
		u-vel:mesh = "adcirc_mesh" ;
		u-vel:units = "m s-1" ;
		u-vel:positive = "east" ;
		u-vel:dry_Value = -99999. ;
	double v-vel(time, node) ;
		v-vel:_FillValue = -99999. ;
		v-vel:long_name = "water column vertically averaged north/south velocity" ;
		v-vel:standard_name = "barotropic_northward_sea_water_velocity" ;
		v-vel:coordinates = "time y x" ;
		v-vel:location = "node" ;
		v-vel:mesh = "adcirc_mesh" ;
		v-vel:units = "m s-1" ;
		v-vel:positive = "north" ;
		v-vel:dry_Value = -99999. ;

// global attributes:
		:primitive_weighting_in_continuity_equation = "unitless" ;
		:sea_surface_height_above_geoid = "meters" ;
		:average_horizontal_eddy_viscosity_in_sea_water_wrt_depth = "(length**2)" ;
}

When loading the files in QGIS/MDAL I got an error:

Following the UGRID Convetions, I wrote a simple python script (this one is an improved v2) to join the zeta and uv files into one output file while fixing some attributes:

netcdf ADCIRC_BG_20200312 {
dimensions:
	time = UNLIMITED ; // (191 currently)
	node = 12769 ;
	nele = 23860 ;
	nvertex = 3 ;
	single = 1 ;
variables:
	double time(time) ;
		time:long_name = "time" ;
		time:standard_name = "time" ;
		time:units = "seconds since 1970-01-01T00:00:00+00:00" ;
		time:calendar = "proleptic_gregorian" ;
	float latitude(node) ;
		latitude:long_name = "latitude" ;
		latitude:standard_name = "latitude" ;
		latitude:units = "degrees_north" ;
	float longitude(node) ;
		longitude:long_name = "longitude" ;
		longitude:standard_name = "longitude" ;
		longitude:units = "degrees_east" ;
	int element(nele, nvertex) ;
		element:long_name = "element" ;
		element:standard_name = "face_node_connectivity" ;
		element:units = "nondimensional" ;
		element:start_index = 1 ;
	int mesh_topology(single) ;
		mesh_topology:long_name = "mesh topology" ;
		mesh_topology:standard_name = "mesh_topology" ;
		mesh_topology:dimension = 2 ;
		mesh_topology:node_coordinates = "longitude latitude" ;
		mesh_topology:face_node_connectivity = "element" ;
		mesh_topology:cf_role = "mesh_topology" ;
		mesh_topology:topology_dimension = 2 ;
	float bathymetry(node) ;
		bathymetry:_FillValue = NaNf ;
		bathymetry:long_name = "bathymetry" ;
		bathymetry:standard_name = "sea_floor_depth" ;
		bathymetry:units = "m" ;
		bathymetry:location = "node" ;
		bathymetry:mesh = "mesh_topology" ;
		bathymetry:coordinates = "latitude longitude" ;
	float ssh(time, node) ;
		ssh:_FillValue = NaNf ;
		ssh:long_name = "sea surface height" ;
		ssh:standard_name = "sea_surface_height_above_geoid" ;
		ssh:units = "m" ;
		ssh:location = "node" ;
		ssh:mesh = "mesh_topology" ;
		ssh:coordinates = "latitude longitude" ;
	float u(time, node) ;
		u:_FillValue = NaNf ;
		u:long_name = "u component of barotropic current" ;
		u:standard_name = "barotropic_eastward_sea_water_velocity" ;
		u:units = "m s-1" ;
		u:location = "node" ;
		u:mesh = "mesh_topology" ;
		u:coordinates = "latitude longitude" ;
	float v(time, node) ;
		v:_FillValue = NaNf ;
		v:long_name = "v component of barotropic current" ;
		v:standard_name = "barotropic_northward_sea_water_velocity" ;
		v:units = "m s-1" ;
		v:location = "node" ;
		v:mesh = "mesh_topology" ;
		v:coordinates = "latitude longitude" ;

// global attributes:
		:title = "ADCIRC" ;
		:institution = "My Institution" ;
		:source = "ADCIRC - ADvanced CIRCulation Model" ;
		:history = "2020-05-06T02:55:42Z: fix_adcirc.py" ;
		:references = "http://www.myinstitution.com/" ;
		:Conventions = "CF-1.7, UGRID-1.0" ;
		:_CoordinateModelRunDate = "2020-03-12T00:00:00Z" ;
}

With this post-processing and the changes done by @PeterPetrik, the file then could be loaded in QGIS 3.12.2 just fine:

Actually, I should have closed this Issue when QGIS 3.12 was released. Closing it now. Thank you very much @PeterPetrik.

[marceloandrioni]

Just to leave a comment here in case someone is using QGIS/MDAL to visualize ADCIRC output files, the files work perfectly with the new Temporal Controller panel in QGIS 3.14 as shown in the image below.

The Controller makes it very easy to select a single instant or to visualize all the times in a loop.

It is also very easy to export the whole data in png files and to create a gif with it.

Note:
The time information can be inserted in the map using:
View >> Decorations >> Title Label >> Insert an Expression >> format_date(@map_start_time, 'yyyy-MM-dd HH:mm:ssZ')
The way to show the legend in the map area is a little hacky but it works. Create an image of the legend (screenshot of the area) and add to the map using:
View >> Decorations >> Image

Thank you.