fix: issue with reading older datx formats

      Documentation and Readme update

README.md

@@ -13,6 +13,7 @@ Superflat_Scripts package requires the following scripts are installed:
 - h5py
 - numpy 
 - scipy
+- matplotlib
 
 It also uses the following packages which should be available with a basic **Conda** installation
 - math
@@ -22,16 +23,16 @@ It also uses the following packages which should be available with a basic **Con
 ####   Installation commands with Conda
    Optionally create a new environment
    
-      conda create -n superflat numpy scipy h5py 
+      conda create -n superflat numpy scipy h5py matplotlib
    
    or install the required packages in an existing environment  
    
       conda activate myenv
-      conda install numpy scipy h5py 
+      conda install numpy scipy h5py matplotlib
 
    Download the Superflat scripts package, either with git command
    
-      git clone ssh://https://gitlab.synchrotron-soleil.fr/OPTIQUE/leaps/superflat_scripts
+      git clone https://gitlab.synchrotron-soleil.fr/OPTIQUE/leaps/superflat_scripts
 
    either by getting and unpacking the zip file available from the Superflat git repository
       
@@ -49,6 +50,13 @@ After activating your environment and switching to Superflat_scripts directory s
 ## Support
 In case you are missing some functionnality or have suggestions for upgrades, please contact Fançois Polack: Francois.polack@synchrotron-soleil.fr
 
+## known issues
+We observed some internal variability in _.datx_ file of different origins. 
+The present implementation is working with _.datx_ files where Surface is stored contiguously without compression.
+Reading will fail if the Surface dataset is stored in a chunked compressed mode. Please report, with the console output, if you face such an issue.
+
+Some _.datx_ files seem to contain more than one surface dataset. **HeightMap.read_Zygofile** will read the first "Surface" dataset pointed in the **/MetaData** table
+
 ## Contributing
 
 This package was started on sample scripts provided by Uwe Flechsig <uwe.flechsig@psi.ch>. It has been extended by Fançois Polack <francois.polack@synchrotron-soleil.fr>,

heightmap.py

 # -*- coding: utf-8 -*-
 
-'''
-This file is part of the Superflat Evaluation Toolkit
-
-Evaluation of surface height data based on  PSD tools
-This file defines the HeightMap class, which is the main storage element of surface height measurements
+## @file heightmap.py
+# This file provides tools for analysing surface height data with statistical and PSD tools.
+#
+# The file defines the HeightMap class, which is the main storage element of surface height measurements.
+#   @author François Polack, Uwe Flechsig
+#   @copyright 2022, Superflat-PCP project  <mailto:superflat-pcp@synchrotron-soleil.fr>
+#   @version 0.1.0
 
-'''
+## @mainpage Superflat Scripts toolkit
+#   This package offers a collection of scripts for processing Optical surface Metrology data.
+#
+#   It aims at providing standardized tools for evaluating the optical quality of a surface in conformity to the
+#   specifications detailed in the SuperFlat-PCP Request for Tender
+#
+#   @section Heightmap
+#   The file heightmap.py defines the HeightMap class which is an data container for evaluating height data recorded from interferometric measurements.
+#
+#   @section processing
+#   The file processing.py defines functions which can be used to detrend data and calculate the PSD with the Tukey as prescribed in the Superflat-PCP Request for Tender
+#
 
 __author__  = "François Polack, Uwe Flechsig"
 __copyright__ = "2022, Superflat-PCP project"
@@ -92,7 +105,7 @@ class HeightMap(object):
         Invalid data are set to NaN, and then masked to be excluded from computation
 
         `HeightMap.rawZ` is kept unmodified throughout processing to avoid reloading the data from the file\n
-        `HeightMap.z` will be redimensionned to the selected  ROI
+        `HeightMap.z` will be re-dimensioned to the selected  ROI
 
         @b CAUTION  It seems that multiplying a `numpy.ma.array` by a scalar on the right is corrupting the mask,
         while multiplying on the left appears to give a correct result
@@ -105,12 +118,30 @@ class HeightMap(object):
         if Path(filename).suffix==".datx":
             h5file=h5.File(filename, 'r')
             # directly open the Surface dataset
-            dsSurf=h5file["/Measurement/Surface"]
+            # dsSurf=h5file["/Measurement/Surface"] Unfortunately we can only do that on the most recent datx files,
+            # the older ones do not have the Measurement group with links
+            metadata=h5file["/MetaData"] # open the MetaData DS to get the link to the Surface DS
+
+            for i in range(0,len(metadata)):
+                if metadata[i]['Destination'].decode('ascii') == 'Surface':
+                   break
+            if i==len(metadata):
+                print("The Surface group location was not found in file", filename)
+                return False
+            print(metadata[i]['Destination'],": ", metadata[i-1]['Destination']) # .decode('ascii'))
+
+
+            dsSurf=h5file[metadata[i-1]['Destination'].decode('ascii')]
             #    scale_factor=dsSurf.attrs['Interferometric Scale Factor'][0]*dsSurf.attrs['Obliquity Factor'][0]*dsSurf.attrs['Wavelength'][0]
             scale_factor= dsSurf.attrs['Z Converter'][0][2][1]*dsSurf.attrs['Z Converter'][0][2][2]*dsSurf.attrs['Z Converter'][0][2][3]
-            print("no data=",dsSurf.attrs['No Data'][0])
+
             # self.z=np.ma.masked_greater_equal(dsSurf, dsSurf.attrs['No Data'][0]).set_fill_value(np.nan)   # this masks the invalid data but not really changing their value to nan
-            self.z=np.copy(dsSurf)
+            print("\nfile", filename)
+            # print(dsSurf)
+            print("shape=", dsSurf.shape, "type", dsSurf.dtype, "compression", dsSurf.compression, "chunks", dsSurf.chunks)
+            print("no data=",dsSurf.attrs['No Data'][0])
+
+            self.z=dsSurf[()]
             self.z[np.greater_equal(self.z, dsSurf.attrs['No Data'][0]) ]=np.nan  #  invalid data set to NaN
             self.rawZ=self.z=scale_factor*np.ma.masked_invalid(self.z)  # and then masked
             # set to scale
@@ -159,7 +190,7 @@ class HeightMap(object):
         """! defines the scale of coordinate axis of the active `HeightMap` arrays `HeightMap.z` and `HeightMap.ze
         @param offset [optional]  the display origin of the active arrays
 
-        HeightMap.set_scale redefines the HeightMap.x and `HeightMap.y` array which are internally used for displaying an height map.\n
+        HeightMap.set_scale redefines the `HeightMap.x` and `HeightMap.y` array which are internally used for displaying an height map.\n
         This function is  called internally after loading new data or changing the ROI.
 
         User should not need to call this function unless a new origin of the display coordinates is wanted
@@ -179,14 +210,14 @@ class HeightMap(object):
 
     def find_best_toroid(self):
         """! Find the best second degree 2D polynomial fitting the `HeightMap.z` data on the ROI
-            @return The tuple (**parameters**, **residuals**`) where\n
-            `parameters` is the tuple  (**angle**, **rx**, **pitch**, **ry**, **roll**, **Z0**) with:\n
-            **angle** being the angle of rotation of the axes of the best fitting toroid with respect to current ccordinates (in radians)\n
-            **rx** the radius of curvature in the rotated x direction (in m^-1)\n
-            **pitch** the tilt angle at point (0,0) in the rotated x direction (in radians) \n
-            **ry** the radius of curvature in the rotated Y direction (in m^-1)\n
-            **roll** the tilt angle at point (0,0) in the rotated Y direction (in radians) \n
-            **z0**  the piston displacement (in m)\n
+            @return The tuple (**parameters**, **residuals**) where \n
+                **parameters** is the tuple  (**angle**, **rx**, **pitch**, **ry**, **roll**, **Z0**) with:
+                    - **angle** being the angle of rotation of the axes of the best fitting toroid with respect to current ccordinates (in radians)
+                    - **rx** the radius of curvature in the rotated x direction @f$(in m^-1)@f$
+                    - **pitch** the tilt angle at point (0,0) in the rotated x direction (in radians)
+                    - **ry** the radius of curvature in the rotated Y direction @f$(in m^-1)@f$
+                    - **roll** the tilt angle at point (0,0) in the rotated Y direction (in radians)
+                    - **z0**  the piston displacement (in m) \n
                 **residuals**  is a `numpy.ma.array` of the difference between the input data `HeightMap.z` with the best fitting polynomial
 
             HeightMap.find_best_toroid  calls internally the  processing.fit2D_order2 function

processing.py

@@ -25,7 +25,8 @@ def fit2D_poly(xin,yin,z, powers,rot=0):
         @param z: (numpy.ma.array):
             a 2D numpy masked array  constaining the value of the surface to be fitted
         @param powers: (tuplet list)
-            a list of tuplets (xpow, ypow) respectively specifying the  exponent in x and y  of the allowed polynomial terms  \( (x^{x_{pow}} y^{y_{pow}} ) \)
+            a list of tuplets (xpow, ypow) respectively specifying the  exponent in x and y  of the allowed polynomial terms
+            @f$(x^xpow y^ypow)@f$
         @param rot: [optional]
             the rotation angle [in rd] of the coordinate axis on which the polynomal fit is computed
         @return The tuplet (**coeffs**, **residues**)  (list, numpy.ma.array)
@@ -35,6 +36,8 @@ def fit2D_poly(xin,yin,z, powers,rot=0):
         For allowing to fit a surface in the presence of missing or invalid data at some points of the definition grid,
         the function requires a to receive the surface values in a masked array, the mask of which is screening invalid or unwanted values
     """
+
+
     # grid coords
     x, y = np.meshgrid(xin, yin)
     if(rot !=0 ):
@@ -69,7 +72,7 @@ def fit2D_poly(xin,yin,z, powers,rot=0):
 
 def fit_toroid(x,y,z, rot=0):
     """!     Fit a a surface to a second degree 2D polynomial with the following terms only:
-    (1,    x,    y,    x^{2},  y^2) without the xy term.
+    @f$(1,    x,    y,    x^{2},  y^2)@f$ excluding the xy term.
     @param x: (`numpy.array`):
         a 1D array of cooordinates along the fast varying dimension of array z
     @param y: (`numpy.array`):
@@ -78,19 +81,19 @@ def fit_toroid(x,y,z, rot=0):
         a 2D numpy masked array  constaining the value of the surface to be fitted
     @param rot: [optional]
         the rotation angle of the coordinate axis on which the polynomal fit is computed
-    @return (**coeffs** **residues**) : (list, numpy.ma.array`)
+    @return The tuplet (**coeffs**, **residues**) : (list, numpy.ma.array`)
         - **coeffs**: the list of the coefficients of the fitting polynomial, in the order :  ( 1,    x,    y,    x^2,  y^2 )
         - **residues**:  a masked array of the fit residuals, with the same mask as the input `z` array.
 
-    The fonction calls fit2D_poly
+    The fonction calls processing.fit2D_poly
      """
     powers=((0,0),(1,0),(0,1),(2,0),(0,2))
     return fit2D_poly(x,y,z,powers,rot)
 
 def fit2D_order2(x,y,z,rot=0):
     """!
-    Fit a a surface to a 2D polynomial with all second order terms terms
-    $$ 1,    x,    y,    x^2,  xy,  y^2 $$
+    Fits a a surface to a 2D polynomial with all second order terms terms
+    @f$( 1, x, y, x^2, xy, y^2)@f$
     @param x : (`numpy.array`):
         a 1D array of cooordinates along the fast varying dimension of array z
     @param y : (`numpy.array`):
@@ -99,14 +102,11 @@ def fit2D_order2(x,y,z,rot=0):
         a 2D numpy masked array  constaining the value of the surface to be fitted
     @param rot: [optional]
         the rotation angle (in rd] of the coordinate axis on which the polynomal fit is computed
-
-    @return (**coeffs**, **residues**) : (list, numpy.ma.array`)
-
-    - **coeffs**: the list of the coefficients of the fitting polynomial, in the order ( 1,    x,    y,  x^2,  xy,  y^2)
-
+    @return The tuplet (**coeffs**, **residues**) : (list, numpy.ma.array`)
+        - **coeffs**: the list of the coefficients of the fitting polynomial, in the order @f$( 1, x, y, x^2, xy, y^2)@f$
         - **residues**:  a masked array of the fit residuals, with the same mask as the input z array.
 
-    The function calls fit2D_poly
+    The function calls processing.fit2D_poly
     """
 
     powers=((0,0),(1,0),(0,1),(2,0),(1,1),(0,2))
@@ -114,7 +114,7 @@ def fit2D_order2(x,y,z,rot=0):
 
 def fit_toroid1(xin,yin,z):
     """!
-    @deprecated    Same as `fit_toroid`, but without axis rotation and without calling `fit2D_poly`
+    @deprecated    Same as `processing.fit_toroid`, but without axis rotation and without calling `processing.fit2D_poly`
     """
     # grid coords
     x, y = np.meshgrid(xin, yin)