Raster data management — imagery, elevation, and gridded coverage

Raster data represents the world as a continuous grid of cells or pixels, each containing a measured or derived value. In Danish municipal applications, raster data includes aerial photography, satellite imagery, digital elevation models, temperature surfaces, noise maps, and countless other gridded datasets. Understanding how to effectively manage these diverse raster types in GeoPackage enables comprehensive spatial analysis and decision-making.

Types of Raster Data

Imagery Data

Imagery data captures visual information about the Earth's surface through photographs or sensors.

Aerial Photography (Orthophotos)

Aerial photographs corrected for camera tilt and terrain relief, providing a geometrically accurate "map-like" view.

Danish municipal applications:

Planning and development review
Property assessment
Infrastructure inspection
Change detection over time
Public information portals

Characteristics:

True color (RGB) or color infrared (CIR)
Typically 10-25cm ground resolution in Denmark
Large file sizes (municipality coverage = 10-100GB)
JPEG compression suitable
Updated every 1-3 years

Example: Orthophoto specification

Resolution: 10cm (0.1m per pixel)
Spectral bands: Red, Green, Blue (RGB)
Bit depth: 8 bits per band
Format in GeoPackage: JPEG tiles
Coverage: Entire municipality
File size: ~40GB for medium-sized kommune

Satellite Imagery

Images captured by satellites, offering broader coverage but typically lower resolution than aerial photos.

Use cases:

Regional context and overview
Agricultural monitoring
Environmental assessment
Change detection at larger scales

Characteristics:

Resolution: 0.5m (commercial) to 30m (Landsat)
Multiple spectral bands (visible, infrared, thermal)
Regular revisit times (days to weeks)
Can include analysis products (NDVI, land cover)

Scanned Historical Maps

Digitized historical maps providing temporal context.

Applications:

Historical land use research
Cultural heritage documentation
Legal property boundary research
Public history projects

Considerations:

PNG format preserves detail
May require georeferencing
Variable quality and accuracy
Important metadata: original scale, scan resolution, georeferencing method

Elevation Data

Elevation data represents the height of the Earth's surface.

Digital Terrain Model (DTM)

Bare earth elevation, with vegetation and buildings removed.

Danish municipal uses:

Flood risk modeling
Drainage planning
Cut-and-fill calculations
Sight line analysis
Infrastructure design

Characteristics:

0.4m grid resolution (typical Danish standard)
32-bit floating point values (elevation in meters)
Accuracy: ±5-15cm vertical
Derived from LiDAR point clouds
Updated every 3-5 years

Example storage:

Resolution: 0.4m grid
Data type: 32-bit float
Value range: -2.5m to +180m (typical Danish range)
Tile format: PNG with elevation encoding or TIFF
Compression: Minimal (preserve precision)

Digital Surface Model (DSM)

Elevation of the surface including buildings, vegetation, and other features.

Applications:

3D city modeling
Solar potential analysis
Visual impact assessment
Communication tower planning
Wind modeling

Height Models

Difference between DSM and DTM showing above-ground feature heights.

Uses:

Building height extraction
Tree height mapping
Vegetation monitoring
Urban growth analysis

Slope and Aspect

Derived products showing terrain gradient and orientation.

Applications:

Erosion risk assessment
Road planning
Solar panel placement
Agricultural suitability

Example calculation:

Slope: Rate of elevation change (degrees or percent)
Aspect: Compass direction of slope (0-360 degrees)
Storage: 8-bit integer (0-90 for slope degrees, 0-360 for aspect)

Gridded Coverage Data

Continuous surfaces representing various phenomena.

Environmental Data

Temperature surfaces
Precipitation patterns
Air quality indices
Noise levels
Groundwater levels
Soil types and properties

Infrastructure Coverage

Cellular network coverage
Broadband availability
Emergency service response times
Public transit accessibility
Street lighting coverage

Planning and Analysis Surfaces

Population density
Land value surfaces
Development suitability
Risk assessment scores
Service area coverage

Temporal Data

Change detection surfaces
Time series analysis
Seasonal variations
Growth patterns

Storing Raster Data in GeoPackage

GeoPackage handles raster data through its tile mechanism, but different raster types require different approaches.

Approach 1: Standard RGB Imagery

For true color aerial photos and similar imagery.

-- Create tile pyramid for orthophotos
INSERT INTO gpkg_contents (
    table_name, data_type, identifier, description,
    last_change, min_x, min_y, max_x, max_y, srs_id
) VALUES (
    'orthophotos_2024',
    'tiles',
    'Orthophotos 2024',
    'True color aerial photography, 10cm resolution, RGB, captured April-May 2024',
    datetime('now'),
    550000, 6195000, 575000, 6215000,
    25832
);

-- Define tile matrix set and zoom levels
INSERT INTO gpkg_tile_matrix_set VALUES (
    'orthophotos_2024', 25832,
    550000, 6195000, 575000, 6215000
);

-- Zoom levels from overview to full resolution
INSERT INTO gpkg_tile_matrix VALUES 
    ('orthophotos_2024', 10, 32, 32, 256, 256, 3.052, 3.052),
    ('orthophotos_2024', 14, 512, 512, 256, 256, 0.191, 0.191);

-- Create tile table
CREATE TABLE orthophotos_2024 (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    zoom_level INTEGER NOT NULL,
    tile_column INTEGER NOT NULL,
    tile_row INTEGER NOT NULL,
    tile_data BLOB NOT NULL,
    UNIQUE (zoom_level, tile_column, tile_row)
);

-- Tiles stored as JPEG (good compression for photos)
-- Quality 85, RGB, 256x256 pixels

Approach 2: Elevation Data with Float Values

For elevation models requiring precise values.

Method A: Elevation as Grayscale PNG

Encode elevation values in PNG pixel values for single-band data.

-- Create DTM tile pyramid
INSERT INTO gpkg_contents (
    table_name, data_type, identifier, description,
    last_change, min_x, min_y, max_x, max_y, srs_id
) VALUES (
    'dtm_2024',
    'tiles',
    'Digital Terrain Model 2024',
    'Bare earth elevation model, 0.4m resolution, elevations in meters above DVR90',
    datetime('now'),
    550000, 6195000, 575000, 6215000,
    25832
);

INSERT INTO gpkg_tile_matrix_set VALUES (
    'dtm_2024', 25832,
    550000, 6195000, 575000, 6215000
);

-- Fewer zoom levels needed (elevation changes gradually)
INSERT INTO gpkg_tile_matrix VALUES 
    ('dtm_2024', 12, 128, 128, 256, 256, 0.762, 0.762),
    ('dtm_2024', 14, 512, 512, 256, 256, 0.191, 0.191);

CREATE TABLE dtm_2024 (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    zoom_level INTEGER NOT NULL,
    tile_column INTEGER NOT NULL,
    tile_row INTEGER NOT NULL,
    tile_data BLOB NOT NULL,  -- PNG with encoded elevation
    UNIQUE (zoom_level, tile_column, tile_row)
);

-- Add metadata about encoding
CREATE TABLE dtm_2024_metadata (
    property TEXT PRIMARY KEY,
    value TEXT
);

INSERT INTO dtm_2024_metadata VALUES
    ('encoding', 'grayscale_16bit'),
    ('min_elevation', '-2.5'),
    ('max_elevation', '180.0'),
    ('scale_factor', '0.01'),  -- pixel_value * 0.01 = elevation in meters
    ('offset', '-2.5'),
    ('nodata_value', '65535'),
    ('vertical_datum', 'DVR90'),
    ('units', 'meters');

Encoding elevation in PNG:

import numpy as np
from PIL import Image

def encode_elevation_tile(elevation_array, min_elev, max_elev):
    """
    Encode elevation values as 16-bit grayscale PNG
    
    elevation_array: 256x256 numpy array of elevation values (meters)
    Returns: PNG image data as bytes
    """
    # Scale elevation to 0-65535 range
    # Reserve 65535 for NODATA
    scaled = ((elevation_array - min_elev) / (max_elev - min_elev) * 65534).astype(np.uint16)
    
    # Handle NODATA
    scaled[np.isnan(elevation_array)] = 65535
    
    # Create PNG
    img = Image.fromarray(scaled, mode='I;16')
    
    # Save to bytes
    import io
    buffer = io.BytesIO()
    img.save(buffer, format='PNG')
    return buffer.getvalue()

def decode_elevation_tile(png_bytes, min_elev, max_elev):
    """
    Decode PNG back to elevation values
    """
    img = Image.open(io.BytesIO(png_bytes))
    pixel_values = np.array(img, dtype=np.uint16)
    
    # Convert back to elevations
    elevations = (pixel_values / 65534.0) * (max_elev - min_elev) + min_elev
    
    # Restore NODATA
    elevations[pixel_values == 65535] = np.nan
    
    return elevations

Method B: External TIFF with GeoPackage Reference

For very high precision requirements, store original TIFF files and reference them in GeoPackage.

-- Table referencing external raster files
CREATE TABLE dtm_tiles_external (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    tile_name TEXT UNIQUE NOT NULL,
    file_path TEXT NOT NULL,
    min_x REAL, min_y REAL, max_x REAL, max_y REAL,
    min_elevation REAL,
    max_elevation REAL,
    capture_date DATE,
    processing_date DATE
);

INSERT INTO dtm_tiles_external VALUES (
    1,
    'DTM_550_6195',
    './elevation_data/dtm_550000_6195000.tif',
    550000, 6195000, 551000, 6196000,
    12.5, 45.8,
    '2024-03-15',
    '2024-04-20'
);

-- This approach keeps full precision but loses single-file portability

Approach 3: Multi-band Imagery

For imagery with multiple spectral bands (e.g., multispectral satellite imagery).

-- Satellite imagery with multiple bands
INSERT INTO gpkg_contents VALUES (
    'sentinel2_2024',
    'tiles',
    'Sentinel-2 Imagery',
    'Sentinel-2 multispectral imagery, 10m resolution, bands: Red, Green, Blue, NIR',
    datetime('now'),
    550000, 6195000, 575000, 6215000,
    25832
);

-- Standard tile pyramid setup
-- ...

CREATE TABLE sentinel2_2024 (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    zoom_level INTEGER NOT NULL,
    tile_column INTEGER NOT NULL,
    tile_row INTEGER NOT NULL,
    tile_data BLOB NOT NULL,  -- PNG with RGB composite or analysis product
    UNIQUE (zoom_level, tile_column, tile_row)
);

-- Can store different band combinations in different tile tables
CREATE TABLE sentinel2_2024_ndvi (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    zoom_level INTEGER NOT NULL,
    tile_column INTEGER NOT NULL,
    tile_row INTEGER NOT NULL,
    tile_data BLOB NOT NULL,  -- NDVI rendered as color ramp
    UNIQUE (zoom_level, tile_column, tile_row)
);

Approach 4: Categorical Raster Data

For land cover, soil types, zoning, etc.

-- Land cover classification
INSERT INTO gpkg_contents VALUES (
    'land_cover_2024',
    'tiles',
    'Land Cover 2024',
    'Land cover classification: Urban, Agriculture, Forest, Water, Wetland',
    datetime('now'),
    550000, 6195000, 575000, 6215000,
    25832
);

CREATE TABLE land_cover_2024 (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    zoom_level INTEGER NOT NULL,
    tile_column INTEGER NOT NULL,
    tile_row INTEGER NOT NULL,
    tile_data BLOB NOT NULL,  -- PNG with indexed colors
    UNIQUE (zoom_level, tile_column, tile_row)
);

-- Classification legend
CREATE TABLE land_cover_classes (
    class_value INTEGER PRIMARY KEY,
    class_name TEXT NOT NULL,
    description_da TEXT,
    color_hex TEXT,
    display_order INTEGER
);

INSERT INTO land_cover_classes VALUES
    (1, 'Urban', 'Byområde', '#FF0000', 1),
    (2, 'Agriculture', 'Landbrugsareal', '#FFFF00', 2),
    (3, 'Forest', 'Skov', '#00FF00', 3),
    (4, 'Water', 'Vand', '#0000FF', 4),
    (5, 'Wetland', 'Vådområde', '#00FFFF', 5);

Approach 5: Hillshade and Visualization Products

Derived visualization products for terrain representation.

-- Hillshade for terrain visualization
INSERT INTO gpkg_contents VALUES (
    'hillshade_2024',
    'tiles',
    'Terrain Hillshade',
    'Analytical hillshade from DTM, illumination from NW, altitude 45°',
    datetime('now'),
    550000, 6195000, 575000, 6215000,
    25832
);

CREATE TABLE hillshade_2024 (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    zoom_level INTEGER NOT NULL,
    tile_column INTEGER NOT NULL,
    tile_row INTEGER NOT NULL,
    tile_data BLOB NOT NULL,  -- Grayscale PNG
    UNIQUE (zoom_level, tile_column, tile_row)
);

-- Hillshade generation parameters
CREATE TABLE hillshade_parameters (
    parameter TEXT PRIMARY KEY,
    value TEXT
);

INSERT INTO hillshade_parameters VALUES
    ('azimuth', '315'),  -- Light from northwest
    ('altitude', '45'),   -- 45 degrees above horizon
    ('z_factor', '1.0'),  -- Vertical exaggeration
    ('algorithm', 'Horn'),
    ('source_dtm', 'dtm_2024');

Raster Data Workflows

Workflow 1: Processing Aerial Photography

Complete workflow from source imagery to GeoPackage tiles.

Step 1: Prepare source imagery

# Check source file information
gdalinfo source_orthophoto.tif

# Output example:
# Size: 250000 x 200000 pixels
# Pixel Size: 0.10m x 0.10m
# Coordinate System: EPSG:25832
# Bands: 3 (RGB)

Step 2: Create GeoPackage and define pyramid

-- Create GeoPackage structure (as shown in previous sections)
-- Define gpkg_contents, gpkg_tile_matrix_set, gpkg_tile_matrix

Step 3: Generate tiles using GDAL

# Convert to GeoPackage tiles
gdal_translate \
    -of GPKG \
    -co TILE_FORMAT=JPEG \
    -co QUALITY=85 \
    -co BLOCKXSIZE=256 \
    -co BLOCKYSIZE=256 \
    -co RASTER_TABLE=orthophotos_2024 \
    source_orthophoto.tif \
    municipality.gpkg

# Build pyramid (overview levels)
gdaladdo \
    -r average \
    -ro \  # Read-only mode
    --config COMPRESS_OVERVIEW JPEG \
    --config JPEG_QUALITY_OVERVIEW 85 \
    municipality.gpkg \
    2 4 8 16 32  # Overview factors

Step 4: Verify and optimize

-- Count tiles per zoom level
SELECT 
    zoom_level,
    COUNT(*) AS tile_count,
    ROUND(SUM(LENGTH(tile_data)) / 1024.0 / 1024.0, 2) AS size_mb
FROM orthophotos_2024
GROUP BY zoom_level;

-- Optimize database
VACUUM;
ANALYZE;

Workflow 2: Creating Elevation Tile Pyramid

Processing DTM data into GeoPackage.

Step 1: Prepare elevation data

# Check DTM file
gdalinfo dtm_source.tif

# Output:
# Size: 62500 x 50000
# Pixel Size: 0.40m x 0.40m
# Data Type: Float32
# Band 1: Min=-2.5, Max=180.3

Step 2: Process elevation to tiles

import numpy as np
from osgeo import gdal
import sqlite3
from PIL import Image
import io

def create_elevation_tiles(source_tif, gpkg_path, table_name, zoom_level):
    """Generate elevation tiles from source DTM"""
    
    # Open source raster
    ds = gdal.Open(source_tif)
    band = ds.GetRasterBand(1)
    
    # Get elevation range
    stats = band.GetStatistics(True, True)
    min_elev, max_elev = stats[0], stats[1]
    
    print(f"Elevation range: {min_elev:.2f}m to {max_elev:.2f}m")
    
    # Calculate tile parameters
    pixel_size = 0.4  # meters
    tile_size = 256  # pixels
    tile_ground_size = tile_size * pixel_size  # 102.4m
    
    # Get geotransform
    gt = ds.GetGeoTransform()
    origin_x, origin_y = gt[0], gt[3]
    
    # Connect to GeoPackage
    conn = sqlite3.connect(gpkg_path)
    cursor = conn.cursor()
    
    # Calculate number of tiles
    width, height = ds.RasterXSize, ds.RasterYSize
    n_tiles_x = int(np.ceil(width / tile_size))
    n_tiles_y = int(np.ceil(height / tile_size))
    
    print(f"Generating {n_tiles_x} x {n_tiles_y} = {n_tiles_x * n_tiles_y} tiles")
    
    # Generate tiles
    for tile_row in range(n_tiles_y):
        for tile_col in range(n_tiles_x):
            # Calculate pixel offsets
            x_off = tile_col * tile_size
            y_off = tile_row * tile_size
            
            # Read elevation data
            elevation = band.ReadAsArray(
                x_off, y_off,
                min(tile_size, width - x_off),
                min(tile_size, height - y_off)
            ).astype(np.float32)
            
            # Pad if necessary
            if elevation.shape != (tile_size, tile_size):
                padded = np.full((tile_size, tile_size), np.nan, dtype=np.float32)
                padded[:elevation.shape[0], :elevation.shape[1]] = elevation
                elevation = padded
            
            # Encode to 16-bit PNG
            scaled = ((elevation - min_elev) / (max_elev - min_elev) * 65534).astype(np.uint16)
            scaled[np.isnan(elevation)] = 65535
            
            img = Image.fromarray(scaled, mode='I;16')
            buffer = io.BytesIO()
            img.save(buffer, format='PNG', compress_level=6)
            tile_data = buffer.getvalue()
            
            # Insert into database
            cursor.execute(f'''
                INSERT INTO {table_name} 
                (zoom_level, tile_column, tile_row, tile_data)
                VALUES (?, ?, ?, ?)
            ''', (zoom_level, tile_col, tile_row, tile_data))
        
        if tile_row % 10 == 0:
            print(f"Progress: {tile_row}/{n_tiles_y} rows")
            conn.commit()
    
    conn.commit()
    conn.close()
    
    print("Tile generation complete")

# Usage
create_elevation_tiles(
    'dtm_source.tif',
    'municipality.gpkg',
    'dtm_2024',
    14  # Zoom level for 0.4m resolution
)

Step 3: Generate overview tiles

def generate_elevation_overviews(gpkg_path, table_name, from_zoom, to_zoom):
    """Generate lower zoom level tiles from higher zoom tiles"""
    
    conn = sqlite3.connect(gpkg_path)
    cursor = conn.cursor()
    
    # Get elevation encoding parameters
    cursor.execute('''
        SELECT value FROM dtm_2024_metadata 
        WHERE property IN ('min_elevation', 'max_elevation')
    ''')
    min_elev, max_elev = [float(row[0]) for row in cursor.fetchall()]
    
    # For each zoom level from high to low
    for zoom in range(from_zoom - 1, to_zoom - 1, -1):
        print(f"Generating zoom level {zoom}")
        
        # Get matrix dimensions
        cursor.execute('''
            SELECT matrix_width, matrix_height
            FROM gpkg_tile_matrix
            WHERE table_name = ? AND zoom_level = ?
        ''', (table_name, zoom))
        
        matrix_width, matrix_height = cursor.fetchone()
        
        # Generate each overview tile
        for tile_row in range(matrix_height):
            for tile_col in range(matrix_width):
                # Fetch 4 source tiles from next zoom level
                src_zoom = zoom + 1
                src_col = tile_col * 2
                src_row = tile_row * 2
                
                cursor.execute(f'''
                    SELECT tile_column, tile_row, tile_data
                    FROM {table_name}
                    WHERE zoom_level = ?
                      AND tile_column IN (?, ?)
                      AND tile_row IN (?, ?)
                ''', (src_zoom, src_col, src_col + 1, src_row, src_row + 1))
                
                source_tiles = cursor.fetchall()
                
                if not source_tiles:
                    continue  # No source data
                
                # Combine and resample
                combined = np.full((512, 512), np.nan, dtype=np.float32)
                
                for s_col, s_row, s_data in source_tiles:
                    img = Image.open(io.BytesIO(s_data))
                    pixels = np.array(img, dtype=np.uint16)
                    
                    # Decode elevation
                    elevations = (pixels / 65534.0) * (max_elev - min_elev) + min_elev
                    elevations[pixels == 65535] = np.nan
                    
                    # Place in combined array
                    x = (s_col - src_col) * 256
                    y = (s_row - src_row) * 256
                    combined[y:y+256, x:x+256] = elevations
                
                # Resample to 256x256 using averaging
                from scipy.ndimage import zoom
                resampled = zoom(combined, 0.5, order=1)  # Bilinear
                
                # Encode to PNG
                scaled = ((resampled - min_elev) / (max_elev - min_elev) * 65534).astype(np.uint16)
                scaled[np.isnan(resampled)] = 65535
                
                img = Image.fromarray(scaled, mode='I;16')
                buffer = io.BytesIO()
                img.save(buffer, format='PNG', compress_level=6)
                tile_data = buffer.getvalue()
                
                # Insert overview tile
                cursor.execute(f'''
                    INSERT INTO {table_name}
                    (zoom_level, tile_column, tile_row, tile_data)
                    VALUES (?, ?, ?, ?)
                ''', (zoom, tile_col, tile_row, tile_data))
        
        conn.commit()
        print(f"Zoom level {zoom} complete")
    
    conn.close()

# Generate overviews from zoom 14 down to zoom 10
generate_elevation_overviews('municipality.gpkg', 'dtm_2024', 14, 10)

Workflow 3: Creating Derived Products

Generating hillshade, slope, and aspect from elevation.

def generate_hillshade_tiles(dtm_table, hillshade_table, gpkg_path, zoom_level,
                              azimuth=315, altitude=45):
    """Generate hillshade tiles from DTM tiles"""
    import numpy as np
    from osgeo import gdal
    
    conn = sqlite3.connect(gpkg_path)
    cursor = conn.cursor()
    
    # Get DTM parameters
    cursor.execute('''
        SELECT value FROM dtm_2024_metadata 
        WHERE property IN ('min_elevation', 'max_elevation', 'scale_factor')
    ''')
    params = {row[0]: float(row[1]) for row in cursor.fetchall()}
    
    # Get all tiles at this zoom level
    cursor.execute(f'''
        SELECT tile_column, tile_row, tile_data
        FROM {dtm_table}
        WHERE zoom_level = ?
        ORDER BY tile_row, tile_column
    ''', (zoom_level,))
    
    tiles = cursor.fetchall()
    
    for tile_col, tile_row, tile_data in tiles:
        # Decode elevation
        img = Image.open(io.BytesIO(tile_data))
        pixels = np.array(img, dtype=np.uint16)
        elevation = (pixels / 65534.0) * (params['max_elevation'] - params['min_elevation']) + params['min_elevation']
        elevation[pixels == 65535] = np.nan
        
        # Calculate hillshade
        # Simple algorithm (better to use GDAL's DEMProcessing)
        azimuth_rad = np.radians(azimuth)
        altitude_rad = np.radians(altitude)
        
        # Calculate gradients (dx, dy)
        dx = np.gradient(elevation, axis=1)
        dy = np.gradient(elevation, axis=0)
        
        # Calculate slope and aspect
        slope = np.arctan(np.sqrt(dx**2 + dy**2))
        aspect = np.arctan2(-dy, dx)
        
        # Calculate hillshade
        hillshade = 255 * (
            (np.cos(altitude_rad) * np.cos(slope)) +
            (np.sin(altitude_rad) * np.sin(slope) * 
             np.cos(azimuth_rad - aspect))
        )
        
        hillshade = np.clip(hillshade, 0, 255).astype(np.uint8)
        hillshade[np.isnan(elevation)] = 0
        
        # Create PNG
        img = Image.fromarray(hillshade, mode='L')
        buffer = io.BytesIO()
        img.save(buffer, format='PNG')
        tile_data_hillshade = buffer.getvalue()
        
        # Insert hillshade tile
        cursor.execute(f'''
            INSERT INTO {hillshade_table}
            (zoom_level, tile_column, tile_row, tile_data)
            VALUES (?, ?, ?, ?)
        ''', (zoom_level, tile_col, tile_row, tile_data_hillshade))
    
    conn.commit()
    conn.close()
    print(f"Hillshade generation complete for zoom {zoom_level}")

# Generate hillshade
generate_hillshade_tiles('dtm_2024', 'hillshade_2024', 'municipality.gpkg', 14)

Workflow 4: Temporal Analysis

Managing time series of raster data.

-- Store multiple time periods
CREATE TABLE orthophotos_2020 (...);
CREATE TABLE orthophotos_2022 (...);
CREATE TABLE orthophotos_2024 (...);

-- Temporal metadata table
CREATE TABLE imagery_temporal_index (
    table_name TEXT PRIMARY KEY,
    capture_date_start DATE,
    capture_date_end DATE,
    season TEXT,
    year INTEGER,
    description TEXT,
    FOREIGN KEY (table_name) REFERENCES gpkg_contents(table_name)
);

INSERT INTO imagery_temporal_index VALUES
    ('orthophotos_2020', '2020-04-15', '2020-05-20', 'Spring', 2020, '10cm RGB orthophotos'),
    ('orthophotos_2022', '2022-05-01', '2022-05-25', 'Spring', 2022, '10cm RGB orthophotos'),
    ('orthophotos_2024', '2024-04-10', '2024-05-15', 'Spring', 2024, '10cm RGB orthophotos');

-- Query available imagery
SELECT 
    table_name,
    year,
    capture_date_start,
    julianday(capture_date_end) - julianday(capture_date_start) AS capture_duration_days
FROM imagery_temporal_index
ORDER BY year;

-- Find imagery closest to a specific date
SELECT table_name, year, capture_date_start,
       ABS(julianday('2023-05-01') - julianday(capture_date_start)) AS days_difference
FROM imagery_temporal_index
ORDER BY days_difference
LIMIT 1;

Analysis and Extraction

Extracting Values at Points

def extract_elevation_at_point(gpkg_path, dtm_table, x, y, zoom_level=14):
    """Extract elevation value at a specific coordinate"""
    
    conn = sqlite3.connect(gpkg_path)
    cursor = conn.cursor()
    
    # Get tile matrix parameters
    cursor.execute('''
        SELECT 
            tms.min_x, tms.max_y,
            tm.pixel_x_size, tm.pixel_y_size,
            tm.tile_width, tm.tile_height
        FROM gpkg_tile_matrix tm
        JOIN gpkg_tile_matrix_set tms ON tm.table_name = tms.table_name
        WHERE tm.table_name = ? AND tm.zoom_level = ?
    ''', (dtm_table, zoom_level))
    
    origin_x, origin_y, pixel_x, pixel_y, tile_w, tile_h = cursor.fetchone()
    
    # Calculate tile coordinates
    tile_col = int((x - origin_x) / (pixel_x * tile_w))
    tile_row = int((origin_y - y) / (pixel_y * tile_h))
    
    # Get the tile
    cursor.execute(f'''
        SELECT tile_data FROM {dtm_table}
        WHERE zoom_level = ? AND tile_column = ? AND tile_row = ?
    ''', (zoom_level, tile_col, tile_row))
    
    result = cursor.fetchone()
    if not result:
        return None
    
    # Decode tile
    cursor.execute('''
      SELECT value FROM dtm_2024_metadata
      WHERE property IN ('min_elevation', 'max_elevation')
    ''')
    min_elev, max_elev = [float(row[0]) for row in cursor.fetchall()]
    img = Image.open(io.BytesIO(result[0]))
    pixels = np.array(img, dtype=np.uint16)
    elevation = (pixels / 65534.0) * (max_elev - min_elev) + min_elev
    
    # Calculate pixel within tile
    tile_origin_x = origin_x + tile_col * tile_w * pixel_x
    tile_origin_y = origin_y - tile_row * tile_h * pixel_y
    
    pixel_col = int((x - tile_origin_x) / pixel_x)
    pixel_row = int((tile_origin_y - y) / pixel_y)
    
    if 0 <= pixel_row < tile_h and 0 <= pixel_col < tile_w:
        elev_value = elevation[pixel_row, pixel_col]
        return elev_value if not np.isnan(elev_value) else None
    
    return None
    Usage
    elevation = extract_elevation_at_point(
      'municipality.gpkg',
      'dtm_2024',
      565000,  # X coordinate
      6205000,  # Y coordinate
      zoom_level=14
    )
    print(f"Elevation: {elevation:.2f}m")

Calculating Statistics

def calculate_elevation_statistics(gpkg_path, dtm_table, polygon_wkt, zoom_level):
    """Calculate elevation statistics within a polygon"""
    from shapely import wkt
    from shapely.geometry import box
    
    conn = sqlite3.connect(gpkg_path)
    cursor = conn.cursor()
    
    # Parse polygon
    polygon = wkt.loads(polygon_wkt)
    bounds = polygon.bounds  # (minx, miny, maxx, maxy)
    
    # Get tile parameters
    cursor.execute('''
        SELECT 
            tms.min_x, tms.max_y,
            tm.pixel_x_size, tm.pixel_y_size,
            tm.tile_width, tm.tile_height
        FROM gpkg_tile_matrix tm
        JOIN gpkg_tile_matrix_set tms ON tm.table_name = tms.table_name
        WHERE tm.table_name = ? AND tm.zoom_level = ?
    ''', (dtm_table, zoom_level))
    
    origin_x, origin_y, pixel_x, pixel_y, tile_w, tile_h = cursor.fetchone()
    
    # Calculate tile range
    min_col = int((bounds[0] - origin_x) / (pixel_x * tile_w))
    max_col = int((bounds[2] - origin_x) / (pixel_x * tile_w))
    min_row = int((origin_y - bounds[3]) / (pixel_y * tile_h))
    max_row = int((origin_y - bounds[1]) / (pixel_y * tile_h))
    
    # Get encoding parameters
    cursor.execute('''
        SELECT value FROM dtm_2024_metadata
        WHERE property IN ('min_elevation', 'max_elevation')
    ''')
    min_elev, max_elev = [float(row[0]) for row in cursor.fetchall()]
    
    # Collect elevation values
    all_elevations = []
    
    for tile_row in range(min_row, max_row + 1):
        for tile_col in range(min_col, max_col + 1):
            # Get tile
            cursor.execute(f'''
                SELECT tile_data FROM {dtm_table}
                WHERE zoom_level = ? AND tile_column = ? AND tile_row = ?
            ''', (zoom_level, tile_col, tile_row))
            
            result = cursor.fetchone()
            if not result:
                continue
            
            # Decode
            img = Image.open(io.BytesIO(result[0]))
            pixels = np.array(img, dtype=np.uint16)
            elevation = (pixels / 65534.0) * (max_elev - min_elev) + min_elev
            
            # Mask by polygon (simplified - assumes polygon covers entire tiles)
            valid_elevations = elevation[pixels != 65535]
            all_elevations.extend(valid_elevations.flatten())
    
    # Calculate statistics
    elevations = np.array(all_elevations)
    
    stats = {
        'min': float(np.min(elevations)),
        'max': float(np.max(elevations)),
        'mean': float(np.mean(elevations)),
        'median': float(np.median(elevations)),
        'std': float(np.std(elevations)),
        'count': len(elevations)
    }
    
    conn.close()
    return stats

# Usage
polygon_wkt = 'POLYGON((565000 6205000, 566000 6205000, 566000 6206000, 565000 6206000, 565000 6205000))'
stats = calculate_elevation_statistics(
    'municipality.gpkg',
    'dtm_2024',
    polygon_wkt,
    14
)
print(f"Elevation statistics: {stats}")

Best Practices for Raster Data Management

1. Choose Appropriate Resolutions

-- Match zoom levels to source data resolution
-- Don't create finer resolution tiles than source data supports

-- Example for 10cm aerial photos:
-- Zoom 16: 0.095m/pixel (matches 10cm source)
-- Zoom 15: 0.19m/pixel (2x resample)
-- Zoom 14: 0.38m/pixel (4x resample)
-- ...
-- Don't create zoom 17 (0.047m/pixel) - beyond source resolution

2. Organize by Theme and Time

-- Clear naming convention
orthophotos_2024_spring
orthophotos_2024_fall
dtm_2024
dsm_2024
hillshade_2024
land_cover_2024

-- Not:
raster1, raster2, tiles_v3

3. Document Thoroughly

-- Comprehensive metadata
INSERT INTO gpkg_contents (
    table_name, data_type, identifier, description, ...
) VALUES (
    'orthophotos_2024',
    'tiles',
    'Aerial Photos 2024',
    'True color RGB aerial photography. Resolution: 10cm. Captured: April 15-May 20, 2024. Coverage: entire municipality. Contractor: GeoDanmark. Coordinate system: EPSG:25832. Accuracy: horizontal ±15cm, vertical N/A. Format: JPEG quality 85. Bands: Red, Green, Blue. Radiometric: 8-bit per band.',
    ...
);

-- Additional metadata tables
CREATE TABLE raster_metadata (
    table_name TEXT NOT NULL,
    property TEXT NOT NULL,
    value TEXT,
    PRIMARY KEY (table_name, property)
);

INSERT INTO raster_metadata VALUES
    ('orthophotos_2024', 'contractor', 'GeoDanmark'),
    ('orthophotos_2024', 'capture_method', 'Aerial photography'),
    ('orthophotos_2024', 'camera_model', 'UltraCamX'),
    ('orthophotos_2024', 'processing_software', 'Agisoft Metashape'),
    ('orthophotos_2024', 'horizontal_accuracy_cm', '15'),
    ('orthophotos_2024', 'sun_angle_min_deg', '35'),
    ('orthophotos_2024', 'cloud_coverage_max_pct', '5');

4. Validate Data Quality

-- Check for completeness
SELECT 
    tm.table_name,
    tm.zoom_level,
    tm.matrix_width * tm.matrix_height AS expected_tiles,
    COUNT(t.id) AS actual_tiles,
    100.0 * COUNT(t.id) / (tm.matrix_width * tm.matrix_height) AS completeness_pct
FROM gpkg_tile_matrix tm
LEFT JOIN orthophotos_2024 t ON tm.zoom_level = t.zoom_level
WHERE tm.table_name = 'orthophotos_2024'
GROUP BY tm.table_name, tm.zoom_level;

-- Check for anomalies
SELECT 
    zoom_level,
    MIN(LENGTH(tile_data)) AS min_size,
    AVG(LENGTH(tile_data)) AS avg_size,
    MAX(LENGTH(tile_data)) AS max_size
FROM orthophotos_2024
GROUP BY zoom_level;

-- Find outliers (suspiciously small or large tiles)
SELECT 
    zoom_level, tile_column, tile_row,
    LENGTH(tile_data) AS size_bytes
FROM orthophotos_2024
WHERE LENGTH(tile_data) < 5000  -- Unusually small
   OR LENGTH(tile_data) > 150000  -- Unusually large
ORDER BY LENGTH(tile_data);

5. Optimize Storage

-- Monitor storage usage
SELECT 
    c.table_name,
    c.identifier,
    COUNT(t.id) AS tile_count,
    ROUND(SUM(LENGTH(t.tile_data)) / 1024.0 / 1024.0 / 1024.0, 2) AS size_gb,
    ROUND(AVG(LENGTH(t.tile_data)) / 1024.0, 1) AS avg_tile_kb
FROM gpkg_contents c
JOIN orthophotos_2024 t ON c.table_name = 'orthophotos_2024'
WHERE c.data_type = 'tiles'
GROUP BY c.table_name;

-- Identify storage-heavy zoom levels
SELECT 
    zoom_level,
    COUNT(*) AS tiles,
    ROUND(SUM(LENGTH(tile_data)) / 1024.0 / 1024.0, 2) AS size_mb,
    ROUND(100.0 * SUM(LENGTH(tile_data)) / 
        (SELECT SUM(LENGTH(tile_data)) FROM orthophotos_2024), 1) AS pct_of_total
FROM orthophotos_2024
GROUP BY zoom_level
ORDER BY zoom_level;

6. Maintain Audit Trail

-- Track data lineage
CREATE TABLE raster_processing_log (
    log_id INTEGER PRIMARY KEY AUTOINCREMENT,
    table_name TEXT NOT NULL,
    processing_date DATETIME DEFAULT CURRENT_TIMESTAMP,
    operation TEXT NOT NULL,
    source_file TEXT,
    parameters TEXT,  -- JSON
    operator TEXT,
    notes TEXT
);

INSERT INTO raster_processing_log (
    table_name, operation, source_file, parameters, operator, notes
) VALUES (
    'orthophotos_2024',
    'tile_generation',
    'vejle_ortho_2024_utm32.tif',
    '{"tile_format": "JPEG", "quality": 85, "tile_size": 256}',
    'jsmith',
    'Initial tile generation from contractor deliverable'
);

Summary

Raster data management in GeoPackage requires understanding different raster types—imagery, elevation, and gridded coverage—and applying appropriate storage strategies for each. Aerial photography uses JPEG-compressed RGB tiles, elevation models require precision encoding (often as 16-bit PNG), and derived products like hillshade provide visualization support. Effective workflows include preprocessing source data, generating multi-resolution tile pyramids, creating derived analytical products, and implementing robust quality assurance. Proper documentation, validation, and optimization ensure raster datasets remain accessible, performant, and useful for Danish municipal planning, analysis, and decision-making.

Knowledge Check Quiz

Test your understanding of raster data management in GeoPackage:

Question 1: What is the primary advantage of storing aerial photography as JPEG tiles rather than PNG?

a) JPEG supports transparency
b) JPEG provides better compression for photographic imagery
c) JPEG is lossless
d) JPEG supports more colors

Question 2: For a Digital Terrain Model (DTM) with elevation values, why might you use 16-bit PNG encoding?

a) To support transparency
b) To preserve precision of floating-point elevation values
c) Because JPEG doesn't work for elevation
d) To reduce file size

Question 3: What is the typical ground resolution of modern Danish aerial orthophotos?

a) 1 meter
b) 50 centimeters
c) 10 centimeters
d) 1 centimeter

Question 4: What is the difference between a Digital Terrain Model (DTM) and a Digital Surface Model (DSM)?

a) DTM includes buildings and vegetation, DSM is bare earth
b) DTM is bare earth, DSM includes buildings and vegetation
c) DTM is lower resolution than DSM
d) There is no difference

Question 5: For elevation data stored as encoded PNG, what value is typically reserved to represent NODATA?

a) 0
b) -9999
c) 65535 (maximum value for 16-bit)
d) NULL

Question 6: What is hillshade used for?

a) Calculating slope
b) Visualizing terrain relief with simulated lighting
c) Measuring elevation
d) Classifying land cover

Question 7: When creating a tile pyramid from a 10cm resolution source image, which statement is correct?

a) Create as many zoom levels as possible for flexibility
b) Only create zoom levels down to approximately 10cm/pixel
c) Always create exactly 18 zoom levels
d) Zoom levels don't matter for imagery

Question 8: For multispectral satellite imagery with multiple bands, how might you store the data in GeoPackage?

a) Store each band in a separate GeoPackage file
b) GeoPackage cannot store multispectral data
c) Store RGB composites or derived products (like NDVI) as PNG tiles
d) Store as text files with references

Question 9: What is a typical JPEG quality setting that balances file size and visual quality for aerial photography?

a) 100 (maximum)
b) 85
c) 50
d) 25

Question 10: Why might you create fewer zoom levels for elevation data compared to aerial photography?

a) Elevation data is less important
b) Elevation changes gradually, so extreme detail at every zoom isn't necessary
c) Elevation data takes more storage
d) GeoPackage limits elevation zoom levels

Question 11: When storing land cover classification data, what tile format is most appropriate?

a) JPEG for compression
b) PNG with indexed colors matching classification values
c) WebP for modern browsers
d) TIFF for precision

Question 12: What does the term "radiometric resolution" refer to in imagery?

a) The spatial resolution (cm per pixel)
b) The number of bits per pixel/band (e.g., 8-bit, 16-bit)
c) The number of spectral bands
d) The geometric accuracy

Question 13: For managing temporal aerial photography (multiple years), what is a good practice?

a) Overwrite old imagery with new imagery
b) Store each year in a separate tile table with clear naming
c) Mix all years in one table
d) Delete old imagery to save space

Question 14: What Python library is commonly used for reading and processing raster data before tiling?

a) Pandas
b) NumPy only
c) GDAL/OGR
d) Matplotlib

Question 15: When encoding elevation values into PNG pixels, what mathematical operation converts pixel values back to elevations?

a) pixel_value * 100
b) (pixel_value / max_pixel_value) * (max_elev - min_elev) + min_elev
c) pixel_value + elevation_offset
d) sqrt(pixel_value)

Answer Key

b) JPEG provides better compression for photographic imagery - JPEG's lossy compression is optimized for photographs, resulting in much smaller files than PNG with acceptable quality.
b) To preserve precision of floating-point elevation values - 16-bit PNG allows encoding elevation values with sufficient precision by scaling to the 0-65535 range.
c) 10 centimeters - Modern Danish aerial orthophotos typically have 10cm (0.1m) ground resolution, though some older datasets may be 15-25cm.
b) DTM is bare earth, DSM includes buildings and vegetation - DTM represents the bare ground surface, while DSM includes all above-ground features.
c) 65535 (maximum value for 16-bit) - The maximum 16-bit value (65535) is typically reserved for NODATA, leaving 0-65534 for actual elevation values.
b) Visualizing terrain relief with simulated lighting - Hillshade creates a grayscale representation of terrain with simulated sun illumination for visualization.
b) Only create zoom levels down to approximately 10cm/pixel - Creating finer resolution tiles than the source data provides no additional real detail.
c) Store RGB composites or derived products (like NDVI) as PNG tiles - Multi-band data is typically processed into RGB composites or analytical products for tile storage.
b) 85 - JPEG quality 85 provides a good balance between visual quality and file size for most aerial photography applications.
b) Elevation changes gradually, so extreme detail at every zoom isn't necessary - Unlike imagery with fine details, elevation surfaces are continuous and smooth, requiring fewer zoom levels.
b) PNG with indexed colors matching classification values - PNG supports indexed color palettes, perfect for categorical raster data with discrete classes.
b) The number of bits per pixel/band (e.g., 8-bit, 16-bit) - Radiometric resolution defines how many distinct values can be recorded (e.g., 8-bit = 256 values, 16-bit = 65536 values).
b) Store each year in a separate tile table with clear naming - Separate tables with naming like orthophotos_2024, orthophotos_2022 make temporal data management clear and organized.
c) GDAL/OGR - GDAL is the standard library for geospatial raster (and vector) data processing in Python and many other languages.
b) (pixel_value / max_pixel_value) * (max_elev - min_elev) + min_elev - This formula scales the pixel value range back to the original elevation range.

Updated on Jan 19, 2026

Foundations of GeoPackage

Vector Features

Raster Data and Tiles

Metadata and Relationships

Raster data management — imagery, elevation, and gridded coverage

Types of Raster Data

Imagery Data

Elevation Data

Gridded Coverage Data

Storing Raster Data in GeoPackage

Approach 1: Standard RGB Imagery

Approach 2: Elevation Data with Float Values

Approach 3: Multi-band Imagery

Approach 4: Categorical Raster Data

Approach 5: Hillshade and Visualization Products

Raster Data Workflows

Workflow 1: Processing Aerial Photography

Workflow 2: Creating Elevation Tile Pyramid

Workflow 3: Creating Derived Products

Workflow 4: Temporal Analysis

Analysis and Extraction

Extracting Values at Points

Calculating Statistics

Best Practices for Raster Data Management

1. Choose Appropriate Resolutions

2. Organize by Theme and Time

3. Document Thoroughly

4. Validate Data Quality

5. Optimize Storage

6. Maintain Audit Trail

Summary

Knowledge Check Quiz

Answer Key