H3 Raster to Hex: Converting Geospatial Rasters to Hexagonal Grids¶

import os
import h3
import rasterio
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Polygon
from sqlalchemy import create_engine, text
import json
from typing import List, Dict, Tuple
import logging

class RasterToH3Processor:
    def __init__(self, postgres_url: str):
        self.engine = create_engine(postgres_url)
        self.setup_logging()

    def setup_logging(self):
        """Setup logging."""
        logging.basicConfig(level=logging.INFO)
        self.logger = logging.getLogger(__name__)

    def create_h3_hexes(self, bounds: Tuple[float, float, float, float], 
                       resolution: int) -> List[Dict]:
        """Create H3 hexagons for given bounds and resolution."""
        min_lat, min_lon, max_lat, max_lon = bounds

        # Get all H3 hexagons in bounds
        hexes = []

        # Calculate step size for sampling
        lat_step = (max_lat - min_lat) / 100
        lon_step = (max_lon - min_lon) / 100

        for lat in np.arange(min_lat, max_lat, lat_step):
            for lon in np.arange(min_lon, max_lon, lon_step):
                h3_index = h3.latlng_to_cell(lat, lon, resolution)

                if h3_index not in [h['h3_index'] for h in hexes]:
                    # Get hexagon boundary
                    hex_boundary = h3.cell_to_boundary(h3_index)

                    # Create polygon
                    polygon = Polygon([(lon, lat) for lat, lon in hex_boundary])

                    # Get center coordinates
                    center_lat, center_lon = h3.cell_to_latlng(h3_index)

                    hexes.append({
                        'h3_index': h3_index,
                        'resolution': resolution,
                        'lat': center_lat,
                        'lon': center_lon,
                        'geom': polygon.wkt
                    })

        return hexes

    def process_raster_to_h3(self, raster_path: str, data_type: str, 
                            resolution: int = 9) -> None:
        """Process raster data and convert to H3 hexagons."""
        self.logger.info(f"Processing raster: {raster_path}")

        # Open raster
        with rasterio.open(raster_path) as src:
            # Get raster bounds
            bounds = src.bounds
            transform = src.transform

            # Read raster data
            data = src.read(1)  # Read first band

            # Create H3 hexagons for bounds
            hexes = self.create_h3_hexes(bounds, resolution)

            # Store hexes in database
            self.store_hexes(hexes)

            # Process raster values for each hex
            raster_values = self.aggregate_raster_to_hexes(data, transform, hexes)

            # Store raster data
            self.store_raster_data(raster_values, data_type, raster_path)

    def aggregate_raster_to_hexes(self, data: np.ndarray, transform, 
                                hexes: List[Dict]) -> List[Dict]:
        """Aggregate raster values to H3 hexagons."""
        values = []

        for hex_info in hexes:
            h3_index = hex_info['h3_index']

            # Get hexagon boundary
            hex_boundary = h3.cell_to_boundary(h3_index)

            # Convert to raster coordinates
            raster_coords = []
            for lat, lon in hex_boundary:
                row, col = rasterio.transform.rowcol(transform, lon, lat)
                if 0 <= row < data.shape[0] and 0 <= col < data.shape[1]:
                    raster_coords.append((row, col))

            if raster_coords:
                # Create mask for hexagon
                mask = self.create_hex_mask(raster_coords, data.shape)

                # Extract values within hexagon
                hex_values = data[mask]

                # Calculate statistics
                if len(hex_values) > 0:
                    stats = {
                        'h3_index': h3_index,
                        'value': np.mean(hex_values),
                        'min_value': np.min(hex_values),
                        'max_value': np.max(hex_values),
                        'std_value': np.std(hex_values),
                        'count': len(hex_values)
                    }
                    values.append(stats)

        return values

    def create_hex_mask(self, coords: List[Tuple[int, int]], 
                       shape: Tuple[int, int]) -> np.ndarray:
        """Create mask for hexagon in raster."""
        mask = np.zeros(shape, dtype=bool)

        # Simple polygon filling (in practice, use more sophisticated method)
        for row, col in coords:
            if 0 <= row < shape[0] and 0 <= col < shape[1]:
                mask[row, col] = True

        return mask

    def store_hexes(self, hexes: List[Dict]) -> None:
        """Store H3 hexagons in database."""
        with self.engine.connect() as conn:
            for hex_info in hexes:
                conn.execute(text("""
                    INSERT INTO h3_hexes (h3_index, resolution, lat, lon, geom)
                    VALUES (:h3_index, :resolution, :lat, :lon, ST_GeomFromText(:geom, 4326))
                    ON CONFLICT (h3_index) DO NOTHING
                """), hex_info)
            conn.commit()

    def store_raster_data(self, values: List[Dict], data_type: str, 
                         source: str) -> None:
        """Store raster data in database."""
        with self.engine.connect() as conn:
            for value_info in values:
                conn.execute(text("""
                    INSERT INTO raster_data (h3_index, resolution, value, data_type, source, metadata)
                    VALUES (:h3_index, :resolution, :value, :data_type, :source, :metadata)
                """), {
                    'h3_index': value_info['h3_index'],
                    'resolution': 9,  # Default resolution
                    'value': value_info['value'],
                    'data_type': data_type,
                    'source': source,
                    'metadata': json.dumps({
                        'min_value': value_info['min_value'],
                        'max_value': value_info['max_value'],
                        'std_value': value_info['std_value'],
                        'count': value_info['count']
                    })
                })
            conn.commit()

    def process_multiple_rasters(self, raster_configs: List[Dict]) -> None:
        """Process multiple raster files."""
        for config in raster_configs:
            try:
                self.process_raster_to_h3(
                    config['path'],
                    config['data_type'],
                    config.get('resolution', 9)
                )
                self.logger.info(f"Processed {config['path']}")
            except Exception as e:
                self.logger.error(f"Error processing {config['path']}: {e}")

# Example usage
if __name__ == "__main__":
    processor = RasterToH3Processor(
        postgres_url=os.getenv("POSTGRES_URL", "postgresql://postgres:postgres@localhost:5432/h3_raster")
    )

    # Process sample rasters
    raster_configs = [
        {
            'path': '/data/elevation.tif',
            'data_type': 'elevation',
            'resolution': 9
        },
        {
            'path': '/data/ndvi.tif',
            'data_type': 'ndvi',
            'resolution': 9
        },
        {
            'path': '/data/temperature.tif',
            'data_type': 'temperature',
            'resolution': 8
        }
    ]

    processor.process_multiple_rasters(raster_configs)

# notebooks/h3_resolution_analysis.ipynb
import h3
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import seaborn as sns
from sqlalchemy import create_engine

def analyze_h3_resolutions(postgres_url: str):
    """Analyze H3 data across different resolutions."""
    engine = create_engine(postgres_url)

    # Get data for each resolution
    resolutions = [6, 7, 8, 9, 10, 11]
    results = []

    for res in resolutions:
        query = f"""
        SELECT 
            COUNT(*) as hex_count,
            AVG(value) as avg_value,
            STDDEV(value) as std_value,
            MIN(value) as min_value,
            MAX(value) as max_value
        FROM raster_data rd
        JOIN h3_hexes hh ON rd.h3_index = hh.h3_index
        WHERE hh.resolution = {res}
        """

        df = pd.read_sql(query, engine)
        df['resolution'] = res
        results.append(df)

    # Combine results
    combined_df = pd.concat(results, ignore_index=True)

    # Plot resolution analysis
    fig, axes = plt.subplots(2, 2, figsize=(15, 10))

    # Hex count by resolution
    axes[0, 0].bar(combined_df['resolution'], combined_df['hex_count'])
    axes[0, 0].set_title('Hex Count by Resolution')
    axes[0, 0].set_xlabel('Resolution')
    axes[0, 0].set_ylabel('Hex Count')

    # Average value by resolution
    axes[0, 1].plot(combined_df['resolution'], combined_df['avg_value'])
    axes[0, 1].set_title('Average Value by Resolution')
    axes[0, 1].set_xlabel('Resolution')
    axes[0, 1].set_ylabel('Average Value')

    # Standard deviation by resolution
    axes[1, 0].plot(combined_df['resolution'], combined_df['std_value'])
    axes[1, 0].set_title('Standard Deviation by Resolution')
    axes[1, 0].set_xlabel('Resolution')
    axes[1, 0].set_ylabel('Standard Deviation')

    # Value range by resolution
    axes[1, 1].plot(combined_df['resolution'], combined_df['max_value'] - combined_df['min_value'])
    axes[1, 1].set_title('Value Range by Resolution')
    axes[1, 1].set_xlabel('Resolution')
    axes[1, 1].set_ylabel('Value Range')

    plt.tight_layout()
    plt.show()

    return combined_df

# Usage
results = analyze_h3_resolutions("postgresql://postgres:postgres@localhost:5432/h3_raster")