From NOAA Weather GRIB to PostGIS: Building a Python Workflow for Raster-Based Spatial Analysis

Introduction

Weather and climate datasets are often distributed in specialized raster formats such as GRIB. These formats are optimized for operational forecasting but require careful handling when integrating them into geospatial analysis workflows.

This article demonstrates a Python-based workflow to process NOAA National Weather Service precipitation forecast data, store it as a tiled raster in PostgreSQL/PostGIS, and visualize the results.

The animation above shows the variation across the raster bands contained within the NOAA NDFD precipitation probability dataset. Each frame represents a different raster band, allowing a quick visual assessment of spatial changes in forecast conditions.

This work was inspired by the excellent demonstrations of PostGIS Raster workflows by Paul Ramsey, especially his discussions on raster storage, tiling, and spatial database management. This article is not a reproduction of that workflow; instead, it presents a Python-centered implementation using open-source geospatial libraries for raster processing, database loading, and visualization.

Workflow Overview

The complete workflow follows these steps:

NOAA GRIB Dataset
        ↓
Python / rioxarray / GDAL
        ↓
GeoTIFF + NoData Handling
        ↓
Custom NOAA CRS in PostGIS
        ↓
PostGIS Raster
        ↓
Raster Tiling + Spatial Indexing
        ↓
Natural Earth Populated Places
        ↓
Vector → NOAA Raster CRS
        ↓
Spatial Intersection
        ↓
Raster Value Extraction
        ↓
Location Ranking
        ↓
Folium Visualization

Dataset

For this workflow, I used the NOAA National Weather Service National Digital Forecast Database (NDFD) precipitation probability product.

Dataset details:

• Product: ds.pop12.bin
• Variable: 12-hour Probability of Precipitation
• Format: GRIB
• Projection: NOAA Lambert Conformal Conic
• Coverage: CONUS

The dataset was accessed directly from the NOAA server using GDAL’s virtual file system (/https://tgftp.nws.noaa.gov/SL.us008001/ST.opnl/DF.gr2/DC.ndfd/AR.conus/VP.001-003/ds.pop12.bin/).

The original NOAA Lambert Conformal Conic projection was preserved. Forecasts are made every 6 hours, so we have 5 rasters in multiple bands representing different forecast periods. The animation below provides a visual representation of the spatial variation across those bands.

Probability of Precipitation raster across the available forecast bands

Before loading the raster into PostGIS, the missing areas need to be handled correctly.

Simply replacing missing pixels with a number such as 255 is not enough. Without explicitly defining that value as NoData, GIS software may interpret those pixels as valid raster observations.

Creating the Custom NOAA Spatial Reference

The NOAA raster uses a Lambert Conformal Conic projection based on a spherical Earth model. Because this projection does not directly correspond to a standard EPSG coordinate reference system, I used a custom PostGIS SRID, 99000.
from sqlalchemy import text

with engine.begin() as conn:
 conn.execute(
 text(“””
 INSERT INTO spatial_ref_sys
 (srid, auth_name, auth_srid, srtext, proj4text)
 VALUES
 (:srid, ‘CUSTOM’, :srid, :srtext, :proj4text)
 ON CONFLICT (srid) DO NOTHING
 “””),
 {
 “srid”: 99000,
 “srtext”: noaa_wkt,
 “proj4text”: noaa_proj4
 }
 )

Loading the raster in PostGIS

The processed GeoTIFF was loaded into PostgreSQL/PostGIS Raster using raster2pgsql.

raster2pgsql \
-s 99000 \
-t 256x256 \
-I \
-C \
-M \
pop12_lcc.tif \
weather.pop12

The raster is divided into 256 × 256 tiles and stored as rows in the PostGIS Raster table.

The main options are:

• -s 99000 assigns the custom NOAA spatial reference.
• -t 256x256 divides the raster into manageable tiles.
• -I creates a spatial index.
• -C adds raster constraints.
• -M runs VACUUM ANALYZE after loading.

Tiling is particularly important because spatial queries can identify the relevant raster tile rather than processing the entire raster.

Checking the PostGIS Raster

After loading the raster, I used Python to query PostgreSQL and verify the raster tiles.

Get Sandip Rijal’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

For example, the total number of raster tiles can be checked with:

cursor.execute(“””
 SELECT COUNT(*)
 FROM weather.pop12;
“””)
total_tiles = cursor.fetchone()[0]
print(f"Total tile rows in weather.pop12: {total_tiles}")

Visualizing the Raster Tiles

Once the raster is stored in PostGIS, the individual tile footprints can be extracted for visualization.

Because the raster is stored in the native NOAA Lambert projection, the tile footprints are transformed to EPSG:4326 before being displayed in Folium.

PostGIS Raster Tile
↓
ST_Envelope()
↓
ST_Transform(4326)
↓
GeoJSON
↓
Folium

This provides a useful visual check of how the original weather raster has been divided into spatial tiles inside PostgreSQL.

PostGIS Raster tiles transformed to geographic coordinates and displayed with Folium.

Querying a Single Location

The next step is to move from simply storing and visualizing the raster to actually querying it.

Natural Earth populated-place data are downloaded https://www.naturalearthdata.com/downloads/10m-cultural-vectors/ and stored in PostGIS in EPSG:4326, while the NOAA raster uses the custom Lambert Conformal Conic CRS, 99000.

Display the points over the map of US.

Populated place data point over US

For a single location, the point is transformed from EPSG:4326 to the NOAA raster CRS, matched to the corresponding raster tile, and used to extract the forecast value directly from PostGIS Raster.

query = “””
WITH pt AS (
 SELECT ST_Transform(
 ST_SetSRID(
 ST_MakePoint(-80, 35),
 4326
 ),
 99000
 ) AS pt
)
SELECT
 ST_Value(rast, 1, pt.pt) AS pop
FROM weather.pop12
CROSS JOIN pt
WHERE ST_Intersects(
 pt.pt,
 ST_ConvexHull(rast)
);
“””
result = pd.read_sql(query, engine)

The spatial operation is straightforward:

Point in EPSG:4326
↓
Transform to NOAA CRS (99000)
↓
Identify raster tile
↓
Extract forecast value

The queried point and its corresponding raster footprint can then be visualized:

Single geographic point shown relative to the NOAA raster footprint.

This demonstrates the basic connection between a geographic location and a weather raster stored inside PostGIS.

Extending the Query to US Populated Places

The same approach can be extended from one location to many.

Natural Earth populated-place data are filtered to US locations, and the point geometries are transformed to the NOAA raster CRS during the query.

query = “””
WITH places_99000 AS (
 SELECT
 gid,
 name,
 ST_Transform(geom, 99000) AS geom
 FROM weather.populated_places
 WHERE iso_a2 = ‘US’
)
SELECT
 p.gid,
 p.name,
 ST_Value(r.rast, 4, p.geom) AS pop
FROM places_99000 AS p
JOIN weather.pop12 AS r
 ON ST_Intersects(r.rast, p.geom)
WHERE ST_Value(r.rast, 4, p.geom) IS NOT NULL
ORDER BY pop DESC
LIMIT 10;
“””
results = pd.read_sql(query, engine)

The database performs the complete spatial operation:

US populated places
↓
Transform to NOAA CRS
↓
Find intersecting raster tile
↓
Extract forecast value
↓
Rank locations
↓
Return top locations

The important advantage here is that the complete raster does not need to be loaded into Python for the spatial operation. PostGIS performs the transformation, intersection, and raster-value extraction directly in the database.

The resulting locations and their extracted forecast values can then be converted into geographic features and displayed using Folium.

US top 10 cities with highest probability of rainfall in next 24 hours.

This creates a simple connection between the original NOAA forecast raster, the spatial database, geographic locations, and the final map.

Why Use PostGIS Raster?

A conventional workflow might require the raster and vector data to be loaded into Python before performing the spatial analysis. With PostGIS Raster, the raster and vector data can remain inside the spatial database. This approach becomes increasingly useful when working with large raster datasets, many geographic locations, or repeated spatial queries.

Conclusions

The workflow combines direct NOAA data access, rioxarray and GDAL-based processing, NoData handling, GeoTIFF creation, PostGIS Raster storage, raster tiling, spatial indexing, SQL validation, spatial queries, and Folium visualization.

More importantly, the workflow demonstrates how raster and vector data can be analyzed together directly inside a spatial database.

The core spatial operations are:

ST_Transform()
ST_Intersects()
ST_Value()

Together, they allow a weather raster stored in its native projection to be connected with geographic locations and converted into location-based forecast information.