user defined metadata - spatial (#36)

* user defined metadata - spatial

* discover readme

* read_access for xios3

* gen rec

* f debug

* f debug

* xios3 compatability

* revert import try (old python)

* review actions; sensible naming

---------

Co-authored-by: mo-marqh <[email protected]>

xios_examples/write_spatial_reference/Makefile

@@ -0,0 +1,34 @@
+# Make file for the write demonstartion XIOS programme
+# Targets provided our detailed below...
+#
+# all: (default) Build the write programme
+# clean: Delete all final products and working files
+# run: run the programme
+#
+# Environment Variables expected by this MakeFile:
+#
+# FC: mpif90
+# FCFLAGS: -g & include files for netcdf & xios
+# LD_FLAGS: for xios, netcdf, netcdff, stfc++
+# LD_LIBRARY_PATH: for netCDF & XIOS libs
+# XIOS_BINDIR: The directory for XIOS binary files
+
+.PHONY: all, clean, run
+
+all: write
+
+# fortran compilation
+%.o: %.F90
+	$(FC) $(FCFLAGS) -c $<
+
+# fortran linking
+write: write.o
+	$(FC) -o write.exe write.o $(LDFLAGS) \
+		&& ln -fs $(XIOS_BINDIR)/xios_server.exe .
+
+run:
+	mpiexec -n 1 ./write.exe : -n 1 ./xios_server.exe
+
+# cleanup
+clean:
+	rm -f *.exe *.o *.mod *.MOD *.out *.err

xios_examples/write_spatial_reference/Readme.md

@@ -0,0 +1,3 @@
+Demonstration of user customised metadata from XIOS
+
+This includes spatial coordinate reference system metadata definitions.

xios_examples/write_spatial_reference/axis_check.xml

@@ -0,0 +1,20 @@
+<context>
+
+  <calendar type="Gregorian"/>
+
+  <grid_definition>
+    <grid id="oax_grid">
+      <axis id="lon" />
+      <axis id="lat" />
+      <axis id="alt" />
+    </grid>
+
+  </grid_definition>
+
+  <file_definition type="one_file">
+    <file id="din" name="spatial_data_input" output_freq="1ts" mode="read" enabled=".true.">
+      <field id="odatax" name="original_data" grid_ref="oax_grid" operation="instant" read_access=".true." />
+    </file>
+  </file_definition>
+
+</context>

xios_examples/write_spatial_reference/expected_domain_output.cdl

@@ -0,0 +1,98 @@
+netcdf spatial_data_output {
+dimensions:
+	axis_nbounds = 2 ;
+	lon = 5 ;
+	lat = 5 ;
+	alt = 1 ;
+	time_counter = UNLIMITED ; // (1 currently)
+variables:
+	float lat(lat) ;
+		lat:axis = "Y" ;
+		lat:standard_name = "latitude" ;
+		lat:long_name = "Latitude" ;
+		lat:units = "degrees_north" ;
+	float lon(lon) ;
+		lon:axis = "X" ;
+		lon:standard_name = "longitude" ;
+		lon:long_name = "Longitude" ;
+		lon:units = "degrees_east" ;
+	float alt(alt) ;
+		alt:name = "alt" ;
+		alt:standard_name = "altitude" ;
+		alt:units = "metres" ;
+	double time_instant(time_counter) ;
+		time_instant:standard_name = "time" ;
+		time_instant:long_name = "Time axis" ;
+		time_instant:calendar = "gregorian" ;
+		time_instant:units = "seconds since 2022-02-02 12:00:00" ;
+		time_instant:time_origin = "2022-02-02 12:00:00" ;
+		time_instant:bounds = "time_instant_bounds" ;
+	double time_instant_bounds(time_counter, axis_nbounds) ;
+	double time_counter(time_counter) ;
+		time_counter:axis = "T" ;
+		time_counter:standard_name = "time" ;
+		time_counter:long_name = "Time axis" ;
+		time_counter:calendar = "gregorian" ;
+		time_counter:units = "seconds since 2022-02-02 12:00:00" ;
+		time_counter:time_origin = "2022-02-02 12:00:00" ;
+		time_counter:bounds = "time_counter_bounds" ;
+	double time_counter_bounds(time_counter, axis_nbounds) ;
+	double original_data(time_counter, alt, lat, lon) ;
+		original_data:long_name = "Arbitrary data values" ;
+		original_data:units = "1" ;
+		original_data:online_operation = "instant" ;
+		original_data:interval_operation = "1 h" ;
+		original_data:interval_write = "1 h" ;
+		original_data:cell_methods = "time: point" ;
+		original_data:coordinates = "time_instant" ;
+		original_data:grid_mapping = "wgs84_2d:lat lon egm2008:alt" ;
+	short wgs84_2d ;
+		wgs84_2d:long_name = "WGS84 CRS" ;
+		wgs84_2d:online_operation = "once" ;
+		wgs84_2d:_FillValue = -32767s ;
+		wgs84_2d:missing_value = -32767s ;
+		wgs84_2d:coordinates = "" ;
+		wgs84_2d:crs_wkt = "GEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\", MEMBER[\"World Geodetic System 1984 (Transit)\", ID[\"EPSG\",1166]], MEMBER[\"World Geodetic System 1984 (G730)\", ID[\"EPSG\",1152]], MEMBER[\"World Geodetic System 1984 (G873)\", ID[\"EPSG\",1153]], MEMBER[\"World Geodetic System 1984 (G1150)\", ID[\"EPSG\",1154]], MEMBER[\"World Geodetic System 1984 (G1674)\", ID[\"EPSG\",1155]], MEMBER[\"World Geodetic System 1984 (G1762)\", ID[\"EPSG\",1156]], MEMBER[\"World Geodetic System 1984 (G2139)\", ID[\"EPSG\",1309]], MEMBER[\"World Geodetic System 1984 (G2296)\", ID[\"EPSG\",1383]], ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]],ID[\"EPSG\",7030]], ENSEMBLEACCURACY[2],ID[\"EPSG\",6326]],CS[ellipsoidal,2,ID[\"EPSG\",6422]],AXIS[\"Geodetic latitude (Lat)\",north],AXIS[\"Geodetic longitude (Lon)\",east],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9102]],ID[\"EPSG\",4326]]" ;
+	short egm2008 ;
+		egm2008:online_operation = "once" ;
+		egm2008:_FillValue = -32767s ;
+		egm2008:missing_value = -32767s ;
+		egm2008:coordinates = "" ;
+		egm2008:crs_wkt = "VERTCRS[\"EGM2008 height\",VDATUM[\"EGM2008 geoid\",ID[\"EPSG\",1027]],CS[vertical,1,ID[\"EPSG\",6499]],AXIS[\"Gravity-related height (H)\",up],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]],GEOIDMODEL[\"WGS 84 to EGM2008 height (1)\",ID[\"EPSG\",3858]],GEOIDMODEL[\"WGS 84 to EGM2008 height (2)\",ID[\"EPSG\",3859]],ID[\"EPSG\",3855]]" ;
+
+// global attributes:
+		:Conventions = "CF-1.6" ;
+		:timeStamp = "2024-Nov-07 16:13:45 GMT" ;
+		:uuid = "5fc977f2-af3c-4a77-80df-fe8f0ae155cb" ;
+		:name = "Spatial Metadata demonstration" ;
+		:description = "Spatial metadata for coordinates is controlled." ;
+		:title = "Spatial Metadata demonstration" ;
+data:
+
+ lat = 50, 52, 54, 56, 58 ;
+
+ lon = -4, -2, 0, 2, 4 ;
+
+ alt = 50 ;
+
+ time_instant = 27133200 ;
+
+ time_instant_bounds =
+  27133200, 27133200 ;
+
+ time_counter = 27133200 ;
+
+ time_counter_bounds =
+  27133200, 27133200 ;
+
+ original_data =
+  0, 2, 4, 6, 8,
+  2, 4, 6, 8, 10,
+  4, 6, 8, 10, 12,
+  6, 8, 10, 12, 14,
+  8, 10, 12, 14, 16 ;
+
+ wgs84_2d = _ ;
+
+ egm2008 = _ ;
+}

xios_examples/write_spatial_reference/iodef.xml

@@ -0,0 +1,9 @@
+<?xml version="1.0"?>
+<simulation>
+
+<context id="main" src="main.xml"/>
+<context id="axis_check" src="axis_check.xml"/>
+
+<context id="xios" src="xios.xml"/>
+
+</simulation>

xios_examples/write_spatial_reference/main.xml

@@ -0,0 +1,60 @@
+<context>
+
+  <calendar type="Gregorian"/>
+
+  <axis_definition>
+    <axis id="lon" unit="degrees" standard_name="longitude" />
+    <axis id="lat" unit="degrees" standard_name="latitude" />
+    <axis id="alt" unit="metres" standard_name="altitude" />
+
+  </axis_definition>
+
+
+  <domain_definition>
+
+    <domain id="original_domain" type="rectilinear" />
+    <generate_rectilinear_domain/>
+  </domain_definition>
+
+  <grid_definition>
+
+    <grid id="original_grid">
+      <domain domain_ref="original_domain" />
+      <axis axis_ref="alt" />
+    </grid>
+    <grid id="metadata_grid" >
+      <scalar/>
+    </grid>
+
+  </grid_definition>
+
+  <field_definition prec="8">
+    <field id="odata" name="original_data" grid_ref="original_grid" long_name="Arbitrary data values" unit="1"  />
+  </field_definition>
+
+
+  <file_definition type="one_file">
+    <file id="spatial_data_input" output_freq="1ts" mode="read" enabled=".true.">
+      <field id="odatain" name="original_data" grid_ref="original_grid" operation="instant" read_access=".true." />
+    </file>
+    <file id="spatial_data_output" output_freq="1ts">
+      <field_group operation="instant">
+        <field field_ref="odata" >
+          <variable name="grid_mapping" type="string">wgs84_2d:lat lon egm2008:alt</variable>
+        </field>
+      </field_group>
+      <field_group operation="once">
+        <field id="hcrs" name="wgs84_2d" grid_ref="metadata_grid" prec="2" long_name="WGS84 CRS" default_value="-32767">
+          <variable name="crs_wkt" type="string">GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble", MEMBER["World Geodetic System 1984 (Transit)", ID["EPSG",1166]], MEMBER["World Geodetic System 1984 (G730)", ID["EPSG",1152]], MEMBER["World Geodetic System 1984 (G873)", ID["EPSG",1153]], MEMBER["World Geodetic System 1984 (G1150)", ID["EPSG",1154]], MEMBER["World Geodetic System 1984 (G1674)", ID["EPSG",1155]], MEMBER["World Geodetic System 1984 (G1762)", ID["EPSG",1156]], MEMBER["World Geodetic System 1984 (G2139)", ID["EPSG",1309]], MEMBER["World Geodetic System 1984 (G2296)", ID["EPSG",1383]], ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7030]], ENSEMBLEACCURACY[2],ID["EPSG",6326]],CS[ellipsoidal,2,ID["EPSG",6422]],AXIS["Geodetic latitude (Lat)",north],AXIS["Geodetic longitude (Lon)",east],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",4326]]</variable>
+      </field>
+        <field id="vcrs" name="egm2008" grid_ref="metadata_grid" prec="2" default_value="-32767">
+          <variable name="crs_wkt" type="string">VERTCRS["EGM2008 height",VDATUM["EGM2008 geoid",ID["EPSG",1027]],CS[vertical,1,ID["EPSG",6499]],AXIS["Gravity-related height (H)",up],LENGTHUNIT["metre",1,ID["EPSG",9001]],GEOIDMODEL["WGS 84 to EGM2008 height (1)",ID["EPSG",3858]],GEOIDMODEL["WGS 84 to EGM2008 height (2)",ID["EPSG",3859]],ID["EPSG",3855]]</variable>
+      </field>
+      </field_group>
+      <variable name="name" type="string">Spatial Metadata demonstration</variable>
+      <variable name="description" type="string">Spatial metadata for coordinates is controlled.</variable>
+      <variable name="title" type="string">Spatial Metadata demonstration</variable>
+    </file>
+  </file_definition>
+
+</context>

xios_examples/write_spatial_reference/spatial_data_input.cdl

@@ -0,0 +1,37 @@
+netcdf domain_input {
+dimensions:
+        lon = 5 ;
+        lat = 5 ;
+        alt = 1 ;
+variables:
+        float lon(lon) ;
+                lon:standard_name = "longitude" ;
+                lon:units = "degrees";
+        float lat(lat) ;
+                lat:standard_name = "latitude" ;
+                lat:units = "degrees";
+        float alt(alt) ;
+                alt:standard_name = "altitude" ;
+                alt:units = "metres";                
+        double original_data(alt, lat, lon) ;
+                original_data:long_name = "data values" ;
+                original_data:units = "1";
+
+// global attributes:
+                :title = "Input data for XIOS output." ;
+
+data:
+
+ lat = 50, 52, 54, 56, 58 ;
+
+ lon = -4, -2, 0, 2, 4 ;
+
+ alt = 50 ;
+
+ original_data =  0,  2,  4,  6,  8,
+                  2,  4,  6,  8, 10,
+                  4,  6,  8, 10, 12,
+                  6,  8, 10, 12, 14,
+                  8, 10, 12, 14, 16 ;
+
+}

xios_examples/write_spatial_reference/test_write_metadata.py

@@ -0,0 +1,66 @@
+import copy
+import glob
+import netCDF4
+import numpy as np
+import os
+import subprocess
+import unittest
+
+import xios_examples.shared_testing as xshared
+
+this_path = os.path.realpath(__file__)
+this_dir = os.path.dirname(this_path)
+
+class TestWriteMetadata(xshared._TestCase):
+    test_dir = this_dir
+    transient_inputs = ['spatial_data_input.nc']
+    transient_outputs = ['spatial_data_output.nc']
+    executable = './write.exe'
+    rtol = 5e-03
+
+    @classmethod
+    def make_a_write_test(cls, inf, nc_method='cdl_files',
+                          nclients=1, nservers=1):
+        """
+        this function makes a test case and returns it as a test function,
+        suitable to be dynamically added to a TestCase for running.
+
+        """
+        # always copy for value, don't pass by reference.
+        infcp = copy.copy(inf)
+        # expected by the fortran XIOS program's main.xml
+        inputfile = cls.transient_inputs[0]
+        outputfile = cls.transient_outputs[0]
+        def test_write_metadata(self):
+            # create a netCDF file using nc_method
+            cls.make_netcdf(infcp, inputfile, nc_method=nc_method)
+            cls.run_mpi_xios(nclients=nclients, nservers=nservers)
+
+            # load the result netCDF file
+            runfile = '{}/{}'.format(cls.test_dir, outputfile)
+            assert(os.path.exists(runfile))
+            run_cmd = ['ncdump', runfile]
+            ncd = subprocess.run(run_cmd, check=True, capture_output=True)
+
+            with open(f'{this_dir}/expected_domain_output.cdl') as fin:
+                exptd = fin.read()
+            emsg = ''
+            for eline, aline in zip(exptd.split('\n'),
+                                    ncd.stdout.decode("utf-8").split('\n')):
+                if eline != aline:
+                    # skip timestamp and uuid
+                    # until we work out how not to set these
+                    if not (':timeStamp' in eline or ':uuid' in eline):
+                        emsg += eline + '\n' + aline + '\n\tmismatch\n\n'
+
+            self.assertFalse(emsg, msg=emsg)
+        return test_write_metadata
+
+
+f = f'{this_dir}/spatial_data_input.cdl'
+# unique name for the test
+tname = 'test_{}'.format(os.path.splitext(os.path.basename(f))[0])
+# add the test as an attribute (function) to the test class
+
+setattr(TestWriteMetadata, tname,
+        TestWriteMetadata.make_a_write_test(f))