refactor: move spatial io from pyTMD

DAC/calculate_inverse_barometer.py

@@ -223,7 +223,7 @@ def calculate_inverse_barometer(base_dir, MODEL, YEAR=None, RANGE=None,
     gridtheta = (90.0 - gridlat)*np.pi/180.0
 
     # ellipsoidal parameters of WGS84 ellipsoid
-    wgs84 = pyTMD.datum(ellipsoid='WGS84', units='MKS')
+    wgs84 = pyTMD.spatial.datum(ellipsoid='WGS84', units='MKS')
     # semimajor and semiminor axes of the ellipsoid [m]
     a_axis = wgs84.a_axis
     b_axis = wgs84.b_axis

DAC/interp_DAC_ICESat_GLA12.py

@@ -165,8 +165,8 @@ def interp_DAC_ICESat_GLA12(base_dir, INPUT_FILE,
     unique_indices = unique_indices % len(t)
 
     # parameters for Topex/Poseidon and WGS84 ellipsoids
-    topex = pyTMD.datum(ellipsoid='TOPEX', units='MKS')
-    wgs84 = pyTMD.datum(ellipsoid='WGS84', units='MKS')
+    topex = pyTMD.spatial.datum(ellipsoid='TOPEX', units='MKS')
+    wgs84 = pyTMD.spatial.datum(ellipsoid='WGS84', units='MKS')
     # convert from Topex/Poseidon to WGS84 Ellipsoids
     lat_40HZ,elev_40HZ = pyTMD.spatial.convert_ellipsoid(lat_TPX, elev_TPX,
         topex.a_axis, topex.flat, wgs84.a_axis, wgs84.flat,

DAC/interp_IB_response_ICESat2_ATL06.py

@@ -330,7 +330,7 @@ def interp_IB_response_ICESat2(base_dir, INPUT_FILE, MODEL,
     gridtheta = (90.0 - gridlat)*np.pi/180.0
 
     # ellipsoidal parameters of WGS84 ellipsoid
-    wgs84 = pyTMD.datum(ellipsoid='WGS84', units='MKS')
+    wgs84 = pyTMD.spatial.datum(ellipsoid='WGS84', units='MKS')
     # semimajor and semiminor axes of the ellipsoid [m]
     a_axis = wgs84.a_axis
     b_axis = wgs84.b_axis

DAC/interp_IB_response_ICESat2_ATL07.py

@@ -325,7 +325,7 @@ def interp_IB_response_ICESat2(base_dir, INPUT_FILE, MODEL,
     gridtheta = (90.0 - gridlat)*np.pi/180.0
 
     # ellipsoidal parameters of WGS84 ellipsoid
-    wgs84 = pyTMD.datum(ellipsoid='WGS84', units='MKS')
+    wgs84 = pyTMD.spatial.datum(ellipsoid='WGS84', units='MKS')
     # semimajor and semiminor axes of the ellipsoid [m]
     a_axis = wgs84.a_axis
     b_axis = wgs84.b_axis

DAC/interp_IB_response_ICESat2_ATL11.py

@@ -346,7 +346,7 @@ def interp_IB_response_ICESat2(base_dir, INPUT_FILE, MODEL,
     gridtheta = (90.0 - gridlat)*np.pi/180.0
 
     # ellipsoidal parameters of WGS84 ellipsoid
-    wgs84 = pyTMD.datum(ellipsoid='WGS84', units='MKS')
+    wgs84 = pyTMD.spatial.datum(ellipsoid='WGS84', units='MKS')
     # semimajor and semiminor axes of the ellipsoid [m]
     a_axis = wgs84.a_axis
     b_axis = wgs84.b_axis

DAC/interp_IB_response_ICESat_GLA12.py

@@ -338,8 +338,8 @@ def interp_IB_response_ICESat(base_dir, INPUT_FILE, MODEL,
         epoch=timescale.time._j2000_epoch, standard='UTC')
 
     # parameters for Topex/Poseidon and WGS84 ellipsoids
-    topex = pyTMD.datum(ellipsoid='TOPEX', units='MKS')
-    wgs84 = pyTMD.datum(ellipsoid='WGS84', units='MKS')
+    topex = pyTMD.spatial.datum(ellipsoid='TOPEX', units='MKS')
+    wgs84 = pyTMD.spatial.datum(ellipsoid='WGS84', units='MKS')
     # convert from Topex/Poseidon to WGS84 Ellipsoids
     lat_40HZ,elev_40HZ = pyTMD.spatial.convert_ellipsoid(lat_TPX, elev_TPX,
         topex.a_axis, topex.flat, wgs84.a_axis, wgs84.flat,

DEM/MPI_DEM_ICESat_GLA12.py

@@ -415,8 +415,8 @@ def main():
     fv = fileID['Data_40HZ']['Elevation_Surfaces']['d_elev'].attrs['_FillValue']
 
     # parameters for Topex/Poseidon and WGS84 ellipsoids
-    topex = pyTMD.datum(ellipsoid='TOPEX', units='MKS')
-    wgs84 = pyTMD.datum(ellipsoid='WGS84', units='MKS')
+    topex = pyTMD.spatial.datum(ellipsoid='TOPEX', units='MKS')
+    wgs84 = pyTMD.spatial.datum(ellipsoid='WGS84', units='MKS')
     # convert from Topex/Poseidon to WGS84 Ellipsoids
     lat_40HZ,elev_40HZ = pyTMD.spatial.convert_ellipsoid(lat_TPX, elev_TPX,
         topex.a_axis, topex.flat, wgs84.a_axis, wgs84.flat,

DEM/MPI_interpolate_DEM.py

@@ -457,11 +457,11 @@ def main():
 
     # read input file to extract time, spatial coordinates and data
     if (args.format == 'csv') and (comm.rank == 0):
-        dinput = pyTMD.spatial.from_ascii(args.infile,
+        dinput = gz.spatial.from_ascii(args.infile,
             columns=args.variables, header=args.header,
             verbose=args.verbose)
     elif (args.format == 'netCDF4') and (comm.rank == 0):
-        dinput = pyTMD.spatial.from_netCDF4(args.infile,
+        dinput = gz.spatial.from_netCDF4(args.infile,
             xname=args.variables[2], yname=args.variables[1],
             timename=args.variables[0], varname=args.variables[3],
             verbose=args.verbose)
@@ -470,7 +470,7 @@ def main():
         attrib[args.variables[1]] = dinput['attributes']['y']
         attrib[args.variables[0]] = dinput['attributes']['time']
     elif (args.format == 'HDF5') and (comm.rank == 0):
-        dinput = pyTMD.spatial.from_HDF5(args.infile,
+        dinput = gz.spatial.from_HDF5(args.infile,
             xname=args.variables[2], yname=args.variables[1],
             timename=args.variables[0], varname=args.variables[3],
             verbose=args.verbose)
@@ -479,7 +479,7 @@ def main():
         attrib[args.variables[1]] = dinput['attributes']['y']
         attrib[args.variables[0]] = dinput['attributes']['time']
     elif (args.format in ('GTiff','cog')) and (comm.rank == 0):
-        dinput = pyTMD.spatial.from_geotiff(args.infile,
+        dinput = gz.spatial.from_geotiff(args.infile,
             verbose=args.verbose)
         # copy global geotiff attributes for projection and grid parameters
         for att_name in ['projection','wkt','spacing','extent']:
@@ -703,18 +703,18 @@ def main():
 
         # output interpolated DEM to file
         if (args.format == 'csv'):
-            pyTMD.spatial.to_ascii(output, attrib, args.outfile,
+            gz.spatial.to_ascii(output, attrib, args.outfile,
                 delimiter=',', columns=[args.variables[0],
                 args.variables[1],args.variables[2],'dem_h'],
                 verbose=args.verbose)
         elif (args.format == 'netCDF4'):
-            pyTMD.spatial.to_netCDF4(output, attrib, args.outfile,
+            gz.spatial.to_netCDF4(output, attrib, args.outfile,
                 verbose=args.verbose)
         elif (args.format == 'HDF5'):
-            pyTMD.spatial.to_HDF5(output, attrib, args.outfile,
+            gz.spatial.to_HDF5(output, attrib, args.outfile,
                 verbose=args.verbose)
         elif args.format in ('GTiff','cog'):
-            pyTMD.spatial.to_geotiff(output, attrib, args.outfile,
+            gz.spatial.to_geotiff(output, attrib, args.outfile,
                 varname='dem_h', driver=args.format,
                 verbose=args.verbose)
         # change the permissions level to MODE

DEM/interp_ATL14_DEM_ICESat_GLA12.py

@@ -135,8 +135,8 @@ def interp_ATL14_DEM_ICESat(INPUT_FILE,
         i_track_40HZ[map_1HZ_40HZ] = i_track_1HZ[k]
 
     # parameters for Topex/Poseidon and WGS84 ellipsoids
-    topex = pyTMD.datum(ellipsoid='TOPEX', units='MKS')
-    wgs84 = pyTMD.datum(ellipsoid='WGS84', units='MKS')
+    topex = pyTMD.spatial.datum(ellipsoid='TOPEX', units='MKS')
+    wgs84 = pyTMD.spatial.datum(ellipsoid='WGS84', units='MKS')
     # convert from Topex/Poseidon to WGS84 Ellipsoids
     lat_40HZ,elev_40HZ = pyTMD.spatial.convert_ellipsoid(lat_TPX, elev_TPX,
         topex.a_axis, topex.flat, wgs84.a_axis, wgs84.flat,

GZ/calculate_GZ_ICESat2_ATL06.py

@@ -197,7 +197,8 @@ def calculate_GZ_ICESat2(base_dir, INPUT_FILE,
     # get output directory from input file
     if OUTPUT_DIRECTORY is None:
         OUTPUT_DIRECTORY = INPUT_FILE.parent
-    OUTPUT_DIRECTORY.mkdir(mode=MODE, parents=True, exist_ok=True)
+    if not OUTPUT_DIRECTORY.exists():
+        OUTPUT_DIRECTORY.mkdir(mode=MODE, parents=True, exist_ok=True)
     # set the hemisphere flag based on ICESat-2 granule
     HEM = set_hemisphere(GRAN)
 

GZ/calculate_GZ_ICESat2_ATL11.py

@@ -191,7 +191,8 @@ def calculate_GZ_ICESat2(base_dir, INPUT_FILE,
     if OUTPUT_DIRECTORY is None:
         OUTPUT_DIRECTORY = INPUT_FILE.parent
     # create output directory if it doesn't exist
-    OUTPUT_DIRECTORY.mkdir(mode=MODE, parents=True, exist_ok=True)
+    if not OUTPUT_DIRECTORY.exists():
+        OUTPUT_DIRECTORY.mkdir(mode=MODE, parents=True, exist_ok=True)
     # file format for associated auxiliary files
     file_format = '{0}_{1}_{2}_{3}{4}_{5}{6}_{7}_{8}{9}.h5'
     # set the hemisphere flag based on ICESat-2 granule

GZ/calculate_GZ_ICESat_GLA12.py

@@ -228,7 +228,8 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
     if OUTPUT_DIRECTORY is None:
         OUTPUT_DIRECTORY = INPUT_FILE.parent
     # create output directory if it doesn't exist
-    OUTPUT_DIRECTORY.mkdir(mode=MODE, parents=True, exist_ok=True)
+    if not OUTPUT_DIRECTORY.exists():
+        OUTPUT_DIRECTORY.mkdir(mode=MODE, parents=True, exist_ok=True)
     # full path to output file
     OUTPUT_FILE = OUTPUT_DIRECTORY.joinpath(FILENAME)
     # plot counter

GZ/calculate_grounding_zone.py

@@ -351,26 +351,26 @@ def calculate_grounding_zone(base_dir, input_file, output_file,
     dem_file = f'{input_file.stem}_{DEM_MODEL}{input_file.suffix}'
     f2 = input_file.with_name(dem_file)
     if (FORMAT == 'csv'):
-        d1 = pyTMD.spatial.from_ascii(input_file, columns=VARIABLES,
+        d1 = gz.spatial.from_ascii(input_file, columns=VARIABLES,
             header=HEADER, verbose=VERBOSE)
         v2 = [VARIABLES[0], VARIABLES[1], VARIABLES[2], 'dem_h']
-        d2 = pyTMD.spatial.from_ascii(f2, columns=v2,
+        d2 = gz.spatial.from_ascii(f2, columns=v2,
             header=0, verbose=VERBOSE)
     elif (FORMAT == 'netCDF4'):
-        d1 = pyTMD.spatial.from_netCDF4(input_file,
+        d1 = gz.spatial.from_netCDF4(input_file,
             xname=VARIABLES[2], yname=VARIABLES[1],
             timename=VARIABLES[0], varname=VARIABLES[3],
             verbose=VERBOSE)
-        d2 = pyTMD.spatial.from_netCDF4(f2,
+        d2 = gz.spatial.from_netCDF4(f2,
             xname=VARIABLES[2], yname=VARIABLES[1],
             timename=VARIABLES[0], varname='dem_h',
             verbose=VERBOSE)
     elif (FORMAT == 'HDF5'):
-        d1 = pyTMD.spatial.from_HDF5(input_file,
+        d1 = gz.spatial.from_HDF5(input_file,
             xname=VARIABLES[2], yname=VARIABLES[1],
             timename=VARIABLES[0], varname=VARIABLES[3],
             verbose=VERBOSE)
-        d2 = pyTMD.spatial.from_HDF5(f2,
+        d2 = gz.spatial.from_HDF5(f2,
             xname=VARIABLES[2], yname=VARIABLES[1],
             timename=VARIABLES[0], varname='dem_h',
             verbose=VERBOSE)

doc/environment.yml

@@ -8,14 +8,17 @@ dependencies:
   - graphviz
   - lxml
   - matplotlib
+  - myst-nb
   - numpy
   - numpydoc
   - pip
   - python>=3.6
   - python-dateutil
   - scipy
   - sphinx
+  - sphinx-argparse>=0.4
+  - sphinxcontrib-bibtex
+  - sphinx-design
   - sphinx_rtd_theme
   - pip:
-     - sphinx-argparse>=0.4
      - ..

doc/make.bat

@@ -10,8 +10,6 @@ if "%SPHINXBUILD%" == "" (
 set SOURCEDIR=source
 set BUILDDIR=build
 
-if "%1" == "" goto help
-
 %SPHINXBUILD% >NUL 2>NUL
 if errorlevel 9009 (
 	echo.
@@ -21,10 +19,12 @@ if errorlevel 9009 (
 	echo.may add the Sphinx directory to PATH.
 	echo.
 	echo.If you don't have Sphinx installed, grab it from
-	echo.http://sphinx-doc.org/
+	echo.https://www.sphinx-doc.org/
 	exit /b 1
 )
 
+if "%1" == "" goto help
+
 %SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
 goto end