Added total points computation to categorical stats.

polygonqueryonraster/app.py

@@ -166,76 +166,6 @@ class GeoTIFFService:
         return polygon_mask, combined_raster_data
 
 
-
-    def rasterize_polygon_old(self, polygon, polygon_crs, directory):
-        """
-        Rasterizes a given polygon over the extent of the GeoTIFF tiles after reprojecting it to the CRS of the tiles.
-        It combines data from overlapping tiles and returns the resulting raster data.
-
-        This method dynamically reprojects the polygon to match the CRS of each GeoTIFF tile, ensuring accurate
-        overlay and analysis.
-
-        Args:
-            polygon (Polygon): The polygon to rasterize, in its original CRS.
-            polygon_crs (str): The CRS of the provided polygon, e.g., 'EPSG:4326'.
-            directory (str): The directory containing the GeoTIFF tiles.
-
-        Returns:
-            tuple: A tuple containing the rasterized version of the polygon and the combined raster data from all relevant tiles.
-            The tuple consists of (polygon_mask, combined_raster_data), where `polygon_mask` is a binary mask of the rasterized polygon,
-            and `combined_raster_data` is the aggregated raster data from the tiles overlapping with the polygon.
-
-        Raises:
-            ValueError: If no valid GeoTIFF files are found or if the polygon is outside the bounds of the GeoTIFF tiles.
-        """
-        polygon_mask = None
-        combined_raster_data = None
-        affine_transform = None
-        first_tile = True
-
-        for geotiff_path in glob.glob(os.path.join(directory, '*.tif')):
-            try:
-                with rasterio.open(geotiff_path) as src:
-                    # Create a transformer from the polygon CRS to the tile CRS
-                    transformer = Transformer.from_crs(polygon_crs, src.crs, always_xy=True)
-                    reprojected_polygon = transform(transformer.transform, polygon)
-
-                    minx, miny, maxx, maxy = reprojected_polygon.bounds
-
-                    # Calculate the intersection window
-                    window = from_bounds(minx, miny, maxx, maxy, src.transform)
-
-                    # Check if window has valid dimensions
-                    if window.width > 0 and window.height > 0:
-                        w_affine = src.window_transform(window)
-
-                        # Read all band data in the window
-                        data = src.read(window=window)
-
-                        # Rasterize the polygon into the same shape as the data
-                        poly_mask = rasterize([(polygon, 1)], out_shape=data.shape[1:], transform=w_affine, fill=0, all_touched=True, dtype='uint8')
-
-                        # Merge the mask with previous masks
-                        if first_tile:
-                            polygon_mask = poly_mask
-                            combined_raster_data = data
-                            affine_transform = w_affine
-                            first_tile = False
-                        else:
-                            polygon_mask = np.maximum(polygon_mask, poly_mask)
-                            # Combine the data from different tiles
-                            combined_raster_data = np.dstack((combined_raster_data, data))
-                            # Affine transform should remain the same for all tiles
-
-            except Exception as e:
-                print(f"Error processing file '{geotiff_path}': {e}")
-
-        if polygon_mask is None or combined_raster_data is None or affine_transform is None:
-            raise ValueError("No valid GeoTIFF files found or polygon is outside the GeoTIFF bounds.")
-
-        return polygon_mask, combined_raster_data
-
-
     def compute_stats_on_raster(self, combined_raster_data, polygon_mask, categorical):
         """
         Computes statistics on the raster data within the polygon mask. Handles both categorical
@@ -267,7 +197,11 @@ class GeoTIFFService:
 
                 # Compute categorical statistics for each band
                 unique, counts = np.unique(masked_band_data, return_counts=True)
-                stats_result[f'band_{band}'] = dict(zip(unique.tolist(), counts.tolist()))
+                total_points = np.count_nonzero(masked_band_data)  # Count of non-zero (non-no-data) points
+                stats_result[f'band_{band}'] = {
+                    'total_points': total_points,
+                    'categories': dict(zip(unique.tolist(), counts.tolist()))
+                }
             else:
                 # Compute standard statistics for non-categorical data
                 masked_band = np.ma.masked_values(masked_band_data, 0)  # Assuming 0 is the nodata value