Fix data correctness and API stability issues

- download_sample: use chunked species-slice copy instead of loading
  entire remote dataset into memory, preventing OOM on large stores
- identify_pixel_value: return NaN for NoData instead of 0.0 to
  distinguish absent data from zero biomass measurements
- ZarrStore.from_url: replace deprecated fsspec.get_mapper() with
  zarr.storage.FsspecStore.from_url() for Zarr v3 compatibility
- _calculate_image_size: warn when dimensions are clamped to service
  limits so users know effective resolution is degraded

Author: mihiarc

gridfia/api.py

@@ -1364,9 +1364,6 @@ def download_sample(
         >>> local_path = api.download_sample("durham_nc", "data/durham.zarr")
         >>> results = api.calculate_metrics(local_path)
         """
-        import shutil
-        import urllib.request
-
         if sample not in self._SAMPLE_METADATA:
             available = list(self._SAMPLE_METADATA.keys())
             raise ValueError(
@@ -1378,26 +1375,38 @@ def download_sample(
 
         logger.info(f"Downloading {dataset_info['name']} to {output_path}")
 
-        # For now, use the cloud loading and copy approach
-        # A more efficient implementation would use streaming download
         store = self.load_from_cloud(sample=sample)
 
-        # Copy to local zarr
+        # Copy to local zarr using chunked writes to avoid loading
+        # the entire dataset into memory at once
         output_path.mkdir(parents=True, exist_ok=True)
 
         import zarr
         local_store = zarr.storage.LocalStore(output_path)
         local_root = zarr.open_group(store=local_store, mode='w')
 
-        # Copy arrays and attributes
+        # Copy root attributes
         local_root.attrs.update(store.attrs)
 
-        # Copy biomass array - Zarr v3 doesn't allow both data and dtype
-        import numpy as np
-        biomass_data = np.array(store.biomass[:], dtype=store.dtype)
-        local_root.create_array('biomass', data=biomass_data, chunks=store.chunks)
+        # Pre-allocate the local biomass array with the same shape/chunks/dtype
+        n_species, height, width = store.shape
+        local_biomass = local_root.create_array(
+            'biomass',
+            shape=(n_species, height, width),
+            chunks=store.chunks,
+            dtype=store.dtype,
+        )
+
+        # Copy biomass data one species slice at a time to bound memory usage
+        # Each slice is (1, height, width) which is manageable even for large rasters
+        chunk_size = store.chunks[0] if store.chunks else 1
+        for start in range(0, n_species, chunk_size):
+            end = min(start + chunk_size, n_species)
+            local_biomass[start:end] = store.biomass[start:end]
+            if show_progress:
+                logger.debug(f"Copied species {start}-{end} of {n_species}")
 
-        # Copy species metadata - convert to numpy arrays for Zarr v3 compatibility
+        # Copy species metadata — these are small 1D arrays, safe to materialize
         codes_array = np.array(store.species_codes, dtype='U10')
         names_array = np.array(store.species_names, dtype='U100')
         local_root.create_array('species_codes', data=codes_array)

gridfia/external/fia_client.py

@@ -606,23 +606,32 @@ def get_species_statistics(self, species_code: str) -> Dict:
             )
 
     def identify_pixel_value(
-        self, 
-        species_code: str, 
-        x: float, 
-        y: float, 
+        self,
+        species_code: str,
+        x: float,
+        y: float,
         spatial_ref: str = "102100"
-    ) -> float:
+    ) -> Optional[float]:
         """
         Get biomass value for a species at a specific coordinate.
-        
-        Args:
-            species_code: FIA species code
-            x: X coordinate
-            y: Y coordinate
-            spatial_ref: Spatial reference system (default: Web Mercator)
-            
-        Returns:
-            Biomass value at the location
+
+        Parameters
+        ----------
+        species_code : str
+            FIA species code (e.g., "0131").
+        x : float
+            X coordinate in the specified spatial reference.
+        y : float
+            Y coordinate in the specified spatial reference.
+        spatial_ref : str, default="102100"
+            Spatial reference system (default: Web Mercator).
+
+        Returns
+        -------
+        Optional[float]
+            Biomass value at the location, np.nan if the pixel is NoData
+            (e.g., non-forest or outside coverage), or None if the API
+            response is missing the value key entirely.
         """
         function_name = self._get_function_name(species_code)
         if not function_name:
@@ -640,19 +649,19 @@ def identify_pixel_value(
                 'rasterFunction': function_name
             })
         }
-        
+
         try:
             response = self._rate_limited_request("GET", f"{self.base_url}/identify", params=params)
             response.raise_for_status()
             result = response.json()
-            
+
             if 'value' in result:
                 value = result['value']
                 if value == 'NoData' or value is None:
-                    return 0.0  # No biomass at this location
+                    return float('nan')
                 return float(value)
             return None
-            
+
         except requests.RequestException as e:
             print_error(f"Failed to identify pixel: {e}")
             raise APIConnectionError(
@@ -804,21 +813,44 @@ def _get_function_name(self, species_code: str) -> Optional[str]:
         return None
 
     def _calculate_image_size(
-        self, 
-        bbox: Tuple[float, float, float, float], 
+        self,
+        bbox: Tuple[float, float, float, float],
         pixel_size: float
     ) -> str:
-        """Calculate image size based on bbox and pixel size."""
+        """
+        Calculate image size based on bbox and pixel size.
+
+        If the requested dimensions exceed service limits, the image is
+        clamped and a warning is emitted so users know the effective
+        resolution is coarser than the requested pixel size.
+        """
         width = int((bbox[2] - bbox[0]) / pixel_size)
         height = int((bbox[3] - bbox[1]) / pixel_size)
-        
+
         # Limit to service maximums
         max_width = 15000
         max_height = 4100
-        
+
+        clamped = False
         if width > max_width:
+            clamped = True
+            effective_x_res = (bbox[2] - bbox[0]) / max_width
             width = max_width
         if height > max_height:
+            clamped = True
+            effective_y_res = (bbox[3] - bbox[1]) / max_height
             height = max_height
-            
+
+        if clamped:
+            effective_res = max(
+                (bbox[2] - bbox[0]) / width,
+                (bbox[3] - bbox[1]) / height,
+            )
+            print_warning(
+                f"Requested area exceeds service limits. "
+                f"Image clamped to {width}x{height} pixels "
+                f"(effective resolution ~{effective_res:.1f}m instead of {pixel_size}m). "
+                f"Consider tiling the request for full resolution."
+            )
+
         return f"{width},{height}" 

gridfia/utils/zarr_utils.py

@@ -291,11 +291,12 @@ def from_url(
         storage_options = storage_options or {}
 
         try:
-            # Create filesystem mapper for the URL
-            fs_map = fsspec.get_mapper(url, **storage_options)
+            # Use Zarr v3 FsspecStore instead of deprecated fsspec.get_mapper()
+            fsspec_store = zarr.storage.FsspecStore.from_url(
+                url, storage_options=storage_options, read_only=True
+            )
 
-            # Open as Zarr group with consolidated metadata for efficiency
-            root = zarr.open_group(fs_map, mode='r')
+            root = zarr.open_group(store=fsspec_store, mode='r')
 
             return cls(root, store=None, path=None)
 

tests/integration/test_real_api.py

@@ -154,9 +154,9 @@ def test_real_identify_pixel_value(self, client):
 
         result = client.identify_pixel_value("0316", x, y)  # Red maple
 
-        # Result should be a float or None/0.0 if no data
+        # Result should be a float (possibly NaN for NoData) or None
         assert result is None or isinstance(result, (int, float))
-        if result is not None and result != 0.0:
+        if result is not None and not np.isnan(result):
             assert result >= 0  # Biomass should be non-negative
 
     def test_real_total_biomass_export(self, client, temp_dir):

tests/unit/test_fia_client.py

@@ -689,7 +689,7 @@ def test_identify_pixel_value_custom_spatial_ref(self):
                 assert params['sr'] == '4326'
 
     def test_identify_pixel_value_no_data(self):
-        """Test handling of NoData pixel values."""
+        """Test handling of NoData pixel values returns NaN."""
         client = BigMapRestClient()
 
         with patch.object(client, '_get_function_name', return_value='SPCD_0131_Abies_balsamea'):
@@ -701,10 +701,10 @@ def test_identify_pixel_value_no_data(self):
 
                 result = client.identify_pixel_value('0131', -11500000, 5500000)
 
-                assert result == 0.0  # NoData should return 0.0
+                assert np.isnan(result)  # NoData should return NaN
 
     def test_identify_pixel_value_none_value(self):
-        """Test handling of None pixel values."""
+        """Test handling of None pixel values returns NaN."""
         client = BigMapRestClient()
 
         with patch.object(client, '_get_function_name', return_value='SPCD_0131_Abies_balsamea'):
@@ -716,7 +716,7 @@ def test_identify_pixel_value_none_value(self):
 
                 result = client.identify_pixel_value('0131', -11500000, 5500000)
 
-                assert result == 0.0  # None should return 0.0
+                assert np.isnan(result)  # None should return NaN
 
     def test_identify_pixel_value_no_value_key(self):
         """Test handling when response has no value key."""
@@ -1204,7 +1204,7 @@ def test_real_identify_pixel_value(self):
         try:
             result = client.identify_pixel_value('0131', x, y)
 
-            # Result could be a float value or None/0.0 if no data at location
+            # Result could be a float (possibly NaN for NoData) or None
             assert result is None or isinstance(result, (int, float))
         except Exception as e:
             # Real API tests may fail due to network/service issues