diff --git a/xios_examples/write_spatial_reference/Makefile b/xios_examples/write_spatial_reference/Makefile new file mode 100644 index 0000000..b5e888c --- /dev/null +++ b/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 diff --git a/xios_examples/write_spatial_reference/Readme.md b/xios_examples/write_spatial_reference/Readme.md new file mode 100644 index 0000000..98f112c --- /dev/null +++ b/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. diff --git a/xios_examples/write_spatial_reference/__init__.py b/xios_examples/write_spatial_reference/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/xios_examples/write_spatial_reference/axis_check.xml b/xios_examples/write_spatial_reference/axis_check.xml new file mode 100644 index 0000000..6cc26e3 --- /dev/null +++ b/xios_examples/write_spatial_reference/axis_check.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + + + + + + + + diff --git a/xios_examples/write_spatial_reference/expected_domain_output.cdl b/xios_examples/write_spatial_reference/expected_domain_output.cdl new file mode 100644 index 0000000..62ea1c3 --- /dev/null +++ b/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 = _ ; +} diff --git a/xios_examples/write_spatial_reference/iodef.xml b/xios_examples/write_spatial_reference/iodef.xml new file mode 100644 index 0000000..37acd4d --- /dev/null +++ b/xios_examples/write_spatial_reference/iodef.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/xios_examples/write_spatial_reference/main.xml b/xios_examples/write_spatial_reference/main.xml new file mode 100644 index 0000000..19da736 --- /dev/null +++ b/xios_examples/write_spatial_reference/main.xml @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + wgs84_2d:lat lon egm2008:alt + + + + + 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]] + + + 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]] + + + Spatial Metadata demonstration + Spatial metadata for coordinates is controlled. + Spatial Metadata demonstration + + + + diff --git a/xios_examples/write_spatial_reference/spatial_data_input.cdl b/xios_examples/write_spatial_reference/spatial_data_input.cdl new file mode 100644 index 0000000..8a33dbc --- /dev/null +++ b/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 ; + +} diff --git a/xios_examples/write_spatial_reference/test_write_metadata.py b/xios_examples/write_spatial_reference/test_write_metadata.py new file mode 100644 index 0000000..6fd7fb6 --- /dev/null +++ b/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)) diff --git a/xios_examples/write_spatial_reference/write.F90 b/xios_examples/write_spatial_reference/write.F90 new file mode 100644 index 0000000..2ef4b57 --- /dev/null +++ b/xios_examples/write_spatial_reference/write.F90 @@ -0,0 +1,136 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2025 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- +!> Programme to just write data with coordinate system definition metadata. +!> +program write + use xios + use mpi + + implicit none + + integer :: comm = -1 + integer :: rank = -1 + integer :: npar = 0 + + call initialise() + call simulate() + call finalise() +contains + + subroutine initialise() + + type(xios_date) :: origin + type(xios_date) :: start + type(xios_duration) :: tstep + integer :: mpi_error + integer :: lenx + integer :: leny + integer :: lenz + double precision, dimension (:), allocatable :: latvals, lonvals, altvals + + ! Arbitrary datetime setup, required for XIOS but unused + origin = xios_date(2022, 2, 2, 12, 0, 0) + start = xios_date(2022, 12, 13, 12, 0, 0) + tstep = xios_hour + + ! Initialise MPI and XIOS + call MPI_INIT(mpi_error) + call xios_initialize('client', return_comm=comm) + + call MPI_Comm_rank(comm, rank, mpi_error) + call MPI_Comm_size(comm, npar, mpi_error) + + ! use the axis_check context to obtain sizing information on all arrays + ! for use in defining the main context interpretation + + call xios_context_initialize('axis_check', comm) + call xios_set_time_origin(origin) + call xios_set_start_date(start) + call xios_set_timestep(tstep) + + call xios_close_context_definition() + + ! fetch sizes of axes from the input file for allocate + call xios_get_axis_attr('lon', n_glo=lenx) + call xios_get_axis_attr('lat', n_glo=leny) + call xios_get_axis_attr('alt', n_glo=lenz) + + allocate ( lonvals(lenx) ) + allocate ( latvals(leny) ) + allocate ( altvals(lenz) ) + + ! fetch coordinate value arrays from the input file + call xios_get_axis_attr('lon', value=lonvals) + call xios_get_axis_attr('lat', value=latvals) + call xios_get_axis_attr('alt', value=altvals) + + ! finalise axis_check context, no longer in use + call xios_set_current_context('axis_check') + call xios_context_finalize() + + ! initialize the main context for interacting with the data. + call xios_context_initialize('main', comm) + + call xios_set_time_origin(origin) + call xios_set_start_date(start) + call xios_set_timestep(tstep) + + ! define the horizontal domain and vertical axis using the input file + call xios_set_domain_attr("original_domain", ni_glo=lenx, ni=lenx, nj_glo=leny, nj=leny, ibegin=0, jbegin=0) + call xios_set_domain_attr("original_domain", lonvalue_1d=lonvals, latvalue_1d=latvals) + call xios_set_axis_attr("alt", n_glo=lenz, n=lenz, begin=0) + call xios_set_axis_attr("alt", value=altvals) + call xios_close_context_definition() + + end subroutine initialise + + subroutine finalise() + + integer :: mpi_error + + ! Finalise all XIOS contexts and MPI + call xios_set_current_context('main') + call xios_context_finalize() + + call xios_finalize() + call MPI_Finalize(mpi_error) + + end subroutine finalise + + subroutine simulate() + + type(xios_date) :: current + integer :: ts + integer :: lenx + integer :: leny + integer :: lenz + + ! Allocatable arrays, size is taken from input file + double precision, dimension (:,:,:), allocatable :: inodata + + ! obtain sizing of the grid for the array allocation + + call xios_get_domain_attr('original_domain', ni_glo=lenx) + call xios_get_domain_attr('original_domain', nj_glo=leny) + call xios_get_axis_attr('alt', n_glo=lenz) + + allocate ( inodata(lenz, leny, lenx) ) + + ! Load data from the input file + call xios_recv_field('odatain', inodata) + + do ts=1, 1 + call xios_update_calendar(ts) + call xios_get_current_date(current) + ! Send (copy) the original data to the output file. + call xios_send_field('odata', inodata) + enddo + + deallocate (inodata) + + end subroutine simulate + +end program write diff --git a/xios_examples/write_spatial_reference/xios.xml b/xios_examples/write_spatial_reference/xios.xml new file mode 100644 index 0000000..9ec1df0 --- /dev/null +++ b/xios_examples/write_spatial_reference/xios.xml @@ -0,0 +1,22 @@ + + + + + performance + + + 1.0 + + + + + true + + 100 + + + true + + + +