Make zarr the default output format across the pipeline

- compute_dem_param() now defaults to format='zarr' with zstd compression,
  chunking, and consolidated metadata; falls back to netCDF if format='netcdf'
- compute_horizon() now defaults to format='zarr'
- get_solar_geom() now defaults to format='zarr'
- Topoclass.compute_dem_param(), compute_horizon(), compute_solar_geometry()
  now use zarr by default and load existing zarr stores first, netCDF fallback
- Topoclass.downscale_climate() now defaults to zarr output
  (zarr_store='downscaled.zarr') when not explicitly configured otherwise
- _compute_svf_horayzon() and _compute_horizon_horayzon() fixed to accept dy
  for correct non-square pixel handling

This makes zarr the default persistence format for all intermediate and
final products where zarr support was implemented.

Author: tintino

TopoPyScale/solar_geom.py

@@ -25,7 +25,8 @@ def get_solar_geom(df_position,
                    num_threads=None,
                    fname='ds_solar.nc',
                    project_ouput=Path('./'),
-                   method_solar='nrel_numpy'):
+                   method_solar='nrel_numpy',
+                   format='zarr'):
     """
     Optimized function to compute solar position using vectorized pvlib operations.
     
@@ -129,11 +130,27 @@ def get_solar_geom(df_position,
         'title': 'Solar geometry computed with TopoPyScale',
         'method': method_solar,
         'processing': 'vectorized_pvlib',
-        'performance_note': 'Optimized vectorized calculation'
+        'performance_note': 'Optimized vectorized calculation',
+        'format': format
     })
 
-    print(f'---> Saving to {project_ouput / fname}')
-    te.to_netcdf(ds, fname=project_ouput / fname)
+    output_path = Path(project_ouput)
+    output_path.mkdir(parents=True, exist_ok=True)
+    base_path = output_path / Path(fname).with_suffix('')
+    file_path = base_path.with_suffix('.zarr' if format.lower() == 'zarr' else Path(fname).suffix)
+
+    print(f'---> Saving to {file_path}')
+    if format.lower() == 'zarr':
+        try:
+            from zarr.codecs import BloscCodec
+        except ImportError:
+            from zarr import Blosc as BloscCodec
+        encoding = {v: {'compressor': BloscCodec(cname='zstd', clevel=3, shuffle='bitshuffle', blocksize=0)}
+                    for v in ds.data_vars}
+        ds = ds.chunk({'time': min(744, ds.time.size), 'point_name': min(50, ds.point_name.size)})
+        ds.to_zarr(str(file_path), mode='w', encoding=encoding, zarr_format=3, consolidated=True)
+    else:
+        te.to_netcdf(ds, fname=file_path)
 
     return ds
 

TopoPyScale/topo_param.py

@@ -73,7 +73,7 @@ def _compute_slope_aspect_topocalc(dem_arr, dx, dy):
     """Compute slope and aspect using topocalc backend."""
     return gradient.gradient_d8(dem_arr, dx, dy)
 
-def _compute_svf_horayzon(dem_arr, dx, azimuth_inc=5):
+def _compute_svf_horayzon(dem_arr, dx, dy, azimuth_inc=5):
     """
     Compute sky view factor using HORAYZON backend with optimized parameters.
     
@@ -99,7 +99,7 @@ def _compute_svf_topocalc(dem_arr, dx):
     """Compute sky view factor using topocalc backend."""
     return viewf.viewf(np.double(dem_arr), dx)[0]
 
-def _compute_horizon_horayzon(dem_arr, dx, azimuth, num_threads=None):
+def _compute_horizon_horayzon(dem_arr, dx, dy, azimuth, num_threads=None):
     """
     Compute horizon angles using HORAYZON backend.
     
@@ -356,65 +356,81 @@ def open_dem(dem_file):
 
 
 
-def compute_dem_param(dem_file, fname='ds_param.nc', project_directory=Path('./'), output_folder='outputs'):
+def compute_dem_param(dem_file, fname='ds_param.nc', project_directory=Path('./'),
+                      output_folder='outputs', format='zarr'):
     """
-    Function to compute and derive DEM parameters: slope, aspect, sky view factor
+    Function to compute and derive DEM parameters: slope, aspect, sky view factor.
+
+    Parameters
+    ----------
+    dem_file : str or Path
+        Path to raster file (geotiff). Must be in a local cartesian coordinate system (e.g. UTM).
+    fname : str
+        Output filename (extension ignored when format='zarr'; .zarr is forced).
+    project_directory : Path
+        Project root directory.
+    output_folder : str
+        Subdirectory for outputs.
+    format : {'zarr', 'netcdf'}
+        Output format. Zarr is the default for faster reads and parallel access.
+
+    Returns
+    -------
+    xarray.Dataset
+        Dataset containing elevation, slope, aspect, aspect_cos, aspect_sin, svf.
+    """
+    pdir = Path(project_directory)
 
-    Args:
-        dem_file (str): path to raster file (geotif). Raster must be in local cartesian coordinate system (e.g. UTM)
+    # Force extension based on format
+    base_path = pdir / output_folder / Path(fname).with_suffix('')
+    zarr_path = base_path.with_suffix('.zarr')
+    nc_path = base_path.with_suffix('.nc')
 
-    Returns:  
-        dataset: x, y, elev, slope, aspect, svf
+    file_ds = zarr_path if format.lower() == 'zarr' else nc_path
 
-    """
-    pdir = project_directory
-    file_ds = pdir / output_folder / fname
-    if file_ds.is_file():
-        print(f'\n---> Dataset {fname} found.')
+    if file_ds.exists():
+        print(f'\n---> Dataset {file_ds.name} found.')
         try:
-            ds = xr.open_dataset(file_ds)
-            # Test if we can actually read the elevation data
-            _ = ds.elevation.values
-        except (RuntimeError, OSError) as e:
-            if "NetCDF: HDF error" in str(e) or "filter returned failure" in str(e):
-                print(f'\n---> Dataset {fname} corrupted (HDF/compression error). Trying h5netcdf backend...')
+            if format.lower() == 'zarr':
+                ds = xr.open_zarr(str(file_ds), consolidated=True)
+                ds = ds.compute() if hasattr(ds.elevation.data, 'compute') else ds
+            else:
+                ds = xr.open_dataset(file_ds)
+                # Also try h5netcdf in case of compression issues
                 try:
+                    _ = ds.elevation.values
+                except (RuntimeError, OSError):
                     ds = xr.open_dataset(file_ds, engine='h5netcdf')
                     _ = ds.elevation.values
-                except Exception:
-                    print(f'\n---> h5netcdf backend failed. Regenerating dataset from DEM...')
-                    if Path(dem_file).is_file():
-                        ds = open_dem(dem_file)
-                    else:
-                        raise ValueError(f'ERROR: Dataset corrupted and no DEM available to regenerate')
+        except Exception as e:
+            print(f'\n---> Dataset {file_ds.name} could not be read: {e}')
+            if Path(dem_file).is_file():
+                print(f'\n---> Regenerating dataset from DEM...')
+                ds = open_dem(dem_file)
             else:
-                raise e
-
+                raise ValueError(f'ERROR: Dataset corrupted/unreadable and no DEM available to regenerate')
     else:
         if Path(dem_file).is_file():
-            print(f'\n---> No {fname} Dataset found. DEM {dem_file} available.')
+            print(f'\n---> No {file_ds.name} dataset found. DEM {dem_file} available.')
             ds = open_dem(dem_file)
-
         else:
             raise ValueError(f'ERROR: No DEM or dataset available')
 
     var_in = list(ds.variables.keys())
     print('\n---> Extracting DEM parameters (slope, aspect, svf)')
     dx = ds.x.diff('x').median().values
     dy = ds.y.diff('y').median().values
-    
-    # Safely access elevation data with error handling
+
+    # Safely access elevation data
     try:
         dem_arr = ds.elevation.values
     except (RuntimeError, OSError) as e:
         if "NetCDF: HDF error" in str(e) or "filter returned failure" in str(e):
             print(f'---> Error reading elevation data: {str(e)}')
             print('---> Attempting to load elevation data with different method...')
             try:
-                # Try loading the data chunk by chunk or using compute()
                 dem_arr = ds.elevation.load().values
             except Exception:
-                # If all else fails, regenerate from DEM file
                 print('---> All methods failed. Regenerating from original DEM...')
                 if Path(dem_file).is_file():
                     ds_new = open_dem(dem_file)
@@ -432,7 +448,7 @@ def compute_dem_param(dem_file, fname='ds_param.nc', project_directory=Path('./'
             slope, aspect = _compute_slope_aspect_horayzon(dem_arr, dx, dy)
         else:
             slope, aspect = _compute_slope_aspect_topocalc(dem_arr, dx, dy)
-        
+
         ds['slope'] = (["y", "x"], np.deg2rad(slope))
         ds['aspect'] = (["y", "x"], np.deg2rad(aspect))
         if 'aspect_cos' not in var_in:
@@ -443,13 +459,14 @@ def compute_dem_param(dem_file, fname='ds_param.nc', project_directory=Path('./'
     if 'svf' not in var_in:
         print('Computing svf ...')
         if HORAYZON_AVAILABLE:
-            svf = _compute_svf_horayzon(dem_arr, dx, azimuth_inc=5)
+            svf = _compute_svf_horayzon(dem_arr, dx, dy, azimuth_inc=5)
         else:
             svf = _compute_svf_topocalc(dem_arr, dx)
         ds['svf'] = (["y", "x"], svf)
 
     ds.attrs = dict(description="DEM input parameters to TopoSub",
-                   author="TopoPyScale, https://github.com/ArcticSnow/TopoPyScale")
+                   author="TopoPyScale, https://github.com/ArcticSnow/TopoPyScale",
+                   format=format)
     ds.x.attrs = {'units': 'm'}
     ds.y.attrs = {'units': 'm'}
     ds.elevation.attrs = {'units': 'm'}
@@ -459,19 +476,26 @@ def compute_dem_param(dem_file, fname='ds_param.nc', project_directory=Path('./'
     ds.aspect_sin.attrs = {'units': 'sinus'}
     ds.svf.attrs = {'units': 'ratio', 'standard_name': 'svf', 'long_name': 'Sky view factor'}
 
-    if file_ds.is_file():
-        te.to_netcdf(ds, fname=pdir / output_folder / 'tmp' / fname)
-        ds = None
-        shutil.move(pdir / output_folder / 'tmp' /fname, file_ds)
-        ds = xr.open_dataset(file_ds)
+    print(f'---> Saving DEM parameters to {format}: {file_ds}')
+    if format.lower() == 'zarr':
+        try:
+            from zarr.codecs import BloscCodec
+        except ImportError:
+            from zarr import Blosc as BloscCodec
+        encoding = {v: {'compressor': BloscCodec(cname='zstd', clevel=3, shuffle='bitshuffle', blocksize=0)}
+                    for v in ds.data_vars}
+        # Use reasonable chunk sizes for spatial data
+        ds = ds.chunk({'y': min(512, ds.y.size), 'x': min(512, ds.x.size)})
+        ds.to_zarr(str(file_ds), mode='w', encoding=encoding, zarr_format=3, consolidated=True)
+        ds = xr.open_zarr(str(file_ds), consolidated=True)
     else:
         te.to_netcdf(ds, fname=file_ds)
 
     return ds
 
 
 def compute_horizon(dem_file, azimuth_inc=30, num_threads=None, fname='da_horizon.nc',
-                    output_directory=Path('./outputs'), format='netcdf'):
+                    output_directory=Path('./outputs'), format='zarr'):
     """
     Function to compute horizon angles using the best available backend (HORAYZON preferred).
 
@@ -500,7 +524,7 @@ def compute_horizon(dem_file, azimuth_inc=30, num_threads=None, fname='da_horizo
 
         # HORAYZON can compute all azimuths at once - much faster!
         horizon_cos_elev = _compute_horizon_horayzon(
-            ds.elevation.values, dx, azimuth, num_threads
+            ds.elevation.values, dx, dy, azimuth, num_threads
         )
 
         # Convert back to horizon elevation angles

TopoPyScale/topoclass.py

@@ -178,16 +178,26 @@ def load_project(self):
         else:
             print(f'-> WARNING: Centroid file {self.config.outputs.file.df_centroids} not found')
 
-        if (self.config.outputs.path / self.config.outputs.file.ds_param).is_file():
-            self.toposub.ds_param = xr.open_dataset(
-                self.config.outputs.path / self.config.outputs.file.ds_param)
+        # Load DEM parameters (zarr preferred, netCDF fallback)
+        ds_param_zarr = (self.config.outputs.path / self.config.outputs.file.ds_param).with_suffix('.zarr')
+        ds_param_nc = self.config.outputs.path / self.config.outputs.file.ds_param
+        if ds_param_zarr.is_dir():
+            self.toposub.ds_param = xr.open_zarr(str(ds_param_zarr), consolidated=True)
+            print(f'---> DEM parameter zarr store {ds_param_zarr.name} exists and loaded')
+        elif ds_param_nc.is_file():
+            self.toposub.ds_param = xr.open_dataset(ds_param_nc)
             print(f'---> DEM parameter file {self.config.outputs.file.ds_param} exists and loaded')
         else:
             print(f'-> WARNING: DEM parameter file {self.config.outputs.file.ds_param} not found')
 
-        if (self.config.outputs.path / self.config.outputs.file.ds_solar).is_file():
-            self.ds_solar = xr.open_dataset(
-                self.config.outputs.path / self.config.outputs.file.ds_solar)
+        # Load solar geometry (zarr preferred, netCDF fallback)
+        ds_solar_zarr = (self.config.outputs.path / self.config.outputs.file.ds_solar).with_suffix('.zarr')
+        ds_solar_nc = self.config.outputs.path / self.config.outputs.file.ds_solar
+        if ds_solar_zarr.is_dir():
+            self.ds_solar = xr.open_zarr(str(ds_solar_zarr), consolidated=True)
+            print(f'---> Solar zarr store {ds_solar_zarr.name} exists and loaded')
+        elif ds_solar_nc.is_file():
+            self.ds_solar = xr.open_dataset(ds_solar_nc)
             print(f'---> Solar file {self.config.outputs.file.ds_solar} exists and loaded')
         else:
             print(f'-> WARNING: Solar file {self.config.outputs.file.ds_solar} not found')
@@ -240,11 +250,12 @@ def plot_clusters_map(self, var='cluster_labels', cmap=plt.cm.hsv, figsize=(14,
         def write_landform(self):
             ts.write_landform(self.dem_path, self.ds_param, self.project_directory)
 
-    def compute_dem_param(self):
+    def compute_dem_param(self, format='zarr'):
         self.toposub.ds_param = tp.compute_dem_param(self.config.dem.filepath,
                                                      fname=self.config.outputs.file.ds_param,
                                                      project_directory=self.config.project.directory,
-                                                     output_folder=self.config.outputs.path)
+                                                     output_folder=self.config.outputs.path,
+                                                     format=format)
 
     def search_optimum_number_of_clusters(self,
                                           cluster_range=np.arange(100, 1000, 200),
@@ -548,8 +559,12 @@ def compute_solar_geometry(self):
             for i, start in enumerate(self.time_splitter.start_list):
                 end = self.time_splitter.end_list[i]
                 fname = self.config.outputs.path / self.time_splitter.ds_solar_flist[i]
+                fname_zarr = fname.with_suffix('.zarr')
 
-                if fname.is_file():
+                if fname_zarr.is_dir():
+                    self.ds_solar = xr.open_zarr(str(fname_zarr), consolidated=True)
+                    print(f'---> Solar zarr store {self.time_splitter.ds_solar_flist[i]} exists and loaded')
+                elif fname.is_file():
                     self.ds_solar = xr.open_dataset(fname, chunks='auto', engine='h5netcdf')
                     print(f'---> Solar file {self.time_splitter.ds_solar_flist[i]} exists and loaded')
                 else:
@@ -561,12 +576,17 @@ def compute_solar_geometry(self):
                                                       self.config.project.parallelization.setting.multicore.CPU_cores,
                                                       self.time_splitter.ds_solar_flist[i],
                                                       self.config.outputs.path,
-                                                      method_solar=self.config.solar_position_method)
+                                                      method_solar=self.config.solar_position_method,
+                                                      format='zarr')
 
         else:
-            fname = self.config.outputs.path / self.config.outputs.file.ds_solar
-            if fname.is_file():
-                self.ds_solar = xr.open_dataset(fname, chunks='auto', engine='h5netcdf')
+            ds_solar_zarr = self.config.outputs.path / Path(self.config.outputs.file.ds_solar).with_suffix('.zarr')
+            ds_solar_nc = self.config.outputs.path / self.config.outputs.file.ds_solar
+            if ds_solar_zarr.is_dir():
+                self.ds_solar = xr.open_zarr(str(ds_solar_zarr), consolidated=True)
+                print(f'---> Solar zarr store {ds_solar_zarr.name} exists and loaded')
+            elif ds_solar_nc.is_file():
+                self.ds_solar = xr.open_dataset(ds_solar_nc, chunks='auto', engine='h5netcdf')
                 print(f'---> Solar file {self.config.outputs.file.ds_solar} exists and loaded')
             else:
                 self.ds_solar = sg.get_solar_geom(self.toposub.df_centroids,
@@ -576,9 +596,11 @@ def compute_solar_geometry(self):
                                                   str(self.config.dem.epsg),
                                                   self.config.project.parallelization.setting.multicore.CPU_cores,
                                                   self.config.outputs.file.ds_solar,
-                                                  self.config.outputs.path)
+                                                  self.config.outputs.path,
+                                                  method_solar=self.config.solar_position_method,
+                                                  format='zarr')
 
-    def compute_horizon(self, format='netcdf'):
+    def compute_horizon(self, format='zarr'):
         """
         Function to compute horizon angle and sample values for list of points
 
@@ -628,6 +650,16 @@ def downscale_climate(self):
         f_pattern = self.config.outputs.file.downscaled_pt
         zarr_store = self.config.outputs.file.zarr_store
 
+        # Default output format is zarr
+        if zarr_store and str(zarr_store).strip():
+            pass  # user explicitly configured zarr store
+        else:
+            if f_pattern and str(f_pattern).strip():
+                print(f'---> Defaulting downscaled output to zarr (ignoring configured netCDF file_pattern {f_pattern})')
+            zarr_store = 'downscaled.zarr'
+            self.config.outputs.file.zarr_store = 'downscaled.zarr'
+            f_pattern = None
+
         # Clean up existing outputs for replacement
         if zarr_store is not None and len(str(zarr_store)) > 0:
             zarr_store_path = downscaled_dir / zarr_store