FIRMS Hotspot Data Processor (Python): Streaming and Streamlining Wildfire Monitoring in QGIS

(in under 20 seconds)

1. Automated Data Retrieval: Downloads fire hotspot data from multiple FIRMS sources for the last 7 days, tailored to the current display extent in QGIS.

2. Time Calculation: Calculates the time difference between the current time and the fire detection time, enabling temporal analysis.

3. Data Consolidation: Combines data from all sources into a single CSV file with a timestamp, ensuring easy access and traceability.

4. QGIS Integration: Loads the processed data into QGIS as a new layer, complete with custom styling based on the time since detection.

5. KML Export: Exports KML files categorized by different time ranges of hotspots and the display extent, facilitating further analysis and sharing.

"""

FIRMS Hotspot Data Processor for QGIS

Downloads, processes, and visualizes FIRMS (Fire Information for Resource Management System)

hotspot data within QGIS

Key features:

1. Downloads fire hotspot data from the multiple FIRMS sources for the last 7 days for the current display extent

2. Calculates the time difference between the current time and the fire detection time

4. Creates a combined CSV file with all the processed data

5. Loads the data into QGIS as a new layer with custom styling based on time since detection

6. Exports KML files

CSV file:

The combined CSV file with all processed data is saved in the "~/QGIS_FIRMS_Data" directory.

The filename includes a timestamp: "merged_fires_7days_processed_YYYYMMDD_HHMMSS.csv"

KML files:

Five KML files for different time ranges of hotspots + one kml for display extent

Note: This script requires an active internet connection to fetch data from FIRMS URLs.

"""

import importlib

import subprocess

import sys

def install_package(package):

    subprocess.check_call([sys.executable, "-m", "pip", "install", package])

required_packages = ['requests', 'pytz']

for package in required_packages:

    try:

        importlib.import_module(package)

    except ImportError:

        print(f"{package} not found. Installing...")

        install_package(package)

import logging

import requests

import csv

from datetime import datetime

import os

from qgis.core import *

from qgis.utils import iface

from PyQt5.QtGui import QColor

import pytz

from qgis.core import (

    QgsProject, QgsVectorLayer, QgsSymbol, QgsCoordinateReferenceSystem, 

    QgsCoordinateTransform, QgsRectangle, QgsFeatureRequest, QgsFeature, 

    QgsGeometry, QgsVectorFileWriter, QgsField

)

from qgis.PyQt.QtCore import QVariant

def print_message(message):

    print(f"[FIRMS Hotspot Processor] {message}")

print_message("Script started. Processing FIRMS hotspot data...")

# Set up logging

logging.basicConfig(level=logging.INFO)

logger = logging.getLogger(__name__)

# Directory to save CSV files

save_dir = os.path.expanduser("~/QGIS_FIRMS_Data")

os.makedirs(save_dir, exist_ok=True)

print_message(f"Data will be saved to: {save_dir}")

# URLs for the CSV files with the provided map key

base_urls = [

    "https://firms.modaps.eosdis.nasa.gov/mapserver/wfs/Canada/a8fec73c297a9a43dcedd6c8f98ff089/?SERVICE=WFS&REQUEST=GetFeature&VERSION=2.0.0&TYPENAME=ms:fires_modis_7days&SRSNAME=urn:ogc:def:crs:EPSG::4326&outputformat=csv",

    "https://firms.modaps.eosdis.nasa.gov/mapserver/wfs/Canada/a8fec73c297a9a43dcedd6c8f98ff089/?SERVICE=WFS&REQUEST=GetFeature&VERSION=2.0.0&TYPENAME=ms:fires_snpp_7days&SRSNAME=urn:ogc:def:crs:EPSG::4326&outputformat=csv",

    "https://firms.modaps.eosdis.nasa.gov/mapserver/wfs/Canada/a8fec73c297a9a43dcedd6c8f98ff089/?SERVICE=WFS&REQUEST=GetFeature&VERSION=2.0.0&TYPENAME=ms:fires_noaa20_7days&SRSNAME=urn:ogc:def:crs:EPSG::4326&outputformat=csv",

    "https://firms.modaps.eosdis.nasa.gov/mapserver/wfs/Canada/a8fec73c297a9a43dcedd6c8f98ff089/?SERVICE=WFS&REQUEST=GetFeature&VERSION=2.0.0&TYPENAME=ms:fires_noaa21_7days&SRSNAME=urn:ogc:def:crs:EPSG::4326&outputformat=csv"

]

# Time zone information

utc_zone = pytz.utc

pacific_zone = pytz.timezone('America/Los_Angeles')

def download_and_process_data(base_url, current_time_utc, extent):

    print_message(f"Downloading full dataset from: {base_url}")

    headers = {

        'Cache-Control': 'no-cache',

        'Pragma': 'no-cache'

    }

    response = requests.get(base_url, headers=headers)

    if response.status_code != 200:

        print_message(f"Failed to fetch data. Status code: {response.status_code}")

        return []

    csv_content = response.content.decode('utf-8')

    rows = list(csv.DictReader(csv_content.splitlines()))

    if not rows:

        print_message("No data available from this source.")

        return []

    print_message(f"Downloaded {len(rows)} rows. Filtering for current extent...")

    filtered_rows = [row for row in rows if extent.contains(QgsPointXY(float(row['longitude']), float(row['latitude'])))]

    print_message(f"Processing {len(filtered_rows)} rows of data...")

    for row in filtered_rows:

        if 'acq_date' in row and 'acq_time' in row:

            acq_datetime_str = f"{row['acq_date']} {row['acq_time']}"

            acq_datetime = datetime.strptime(acq_datetime_str, '%Y-%m-%d %H%M')

            acq_datetime_utc = pytz.utc.localize(acq_datetime)

            acq_datetime_pst = acq_datetime_utc.astimezone(pacific_zone)

            # Format the date as MM/DD/YYYY

            formatted_date = acq_datetime_pst.strftime('%m/%d/%Y')

            # Format the time as HH:MM

            formatted_time = acq_datetime_pst.strftime('%H:%M')

            # Combine date, time, and PST indicator

            row['acq_datetime_PST'] = f"{formatted_date} {formatted_time} PST"

            time_difference = current_time_utc - acq_datetime_utc

            row['time_difference_hours'] = time_difference.total_seconds() / 3600.0

    return filtered_rows

print_message("Retrieving current map canvas extent...")

canvas = iface.mapCanvas()

extent = canvas.extent()

if canvas.mapSettings().destinationCrs().authid() != 'EPSG:4326':

    print_message("Transforming extent to EPSG:4326...")

    transform = QgsCoordinateTransform(canvas.mapSettings().destinationCrs(), QgsCoordinateReferenceSystem("EPSG:4326"), QgsProject.instance())

    extent = transform.transformBoundingBox(extent)

print_message("Getting current time in UTC...")

current_time_utc = datetime.now(utc_zone)

print_message("Fetching and processing data from multiple sources...")

all_combined_rows = []

for i, base_url in enumerate(base_urls, 1):

    print_message(f"Processing source {i} of {len(base_urls)}...")

    all_combined_rows.extend(download_and_process_data(base_url, current_time_utc, extent))

print_message(f"Total hotspots found: {len(all_combined_rows)}")

if not all_combined_rows:

    print_message("\n" + "*" * 80)

    print_message("*     NO HOTSPOTS DETECTED IN THE LAST 7 DAYS WITHIN THE CURRENT QGIS MAP EXTENT     *")

    print_message("*" * 80 + "\n")

else:

    print_message("\n" + "*" * 80)

    print_message(f"*     HOTSPOTS FOUND FOR LAST 7 DAYS IN CURRENT DISPLAY EXTENT - {len(all_combined_rows)} HOTSPOTS     *")

    print_message("*" * 80 + "\n")

    timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")

    combined_csv_filename = os.path.join(save_dir, f"merged_fires_7days_processed_{timestamp}.csv")

    print_message(f"Writing data to CSV: {combined_csv_filename}")

    fieldnames = all_combined_rows[0].keys()

    with open(combined_csv_filename, 'w', newline='') as combined_csv_file:

        writer = csv.DictWriter(combined_csv_file, fieldnames=fieldnames)

        writer.writeheader()

        for row in all_combined_rows:

            writer.writerow(row)

    print_message("Loading data into QGIS...")

    uri = f"file:///{os.path.abspath(combined_csv_filename)}?delimiter=,&xField=longitude&yField=latitude&crs=EPSG:4326"

    layer = QgsVectorLayer(uri, f"FIRMS Hotspots/Fires in the last 7 days ({timestamp})", "delimitedtext")

    if not layer.isValid():

        print_message(f"Failed to load layer: {combined_csv_filename}")

    else:

        print_message("Adding layer to QGIS project...")

        QgsProject.instance().addMapLayer(layer)

        # Add the new field to the layer

        layer.startEditing()

        layer.addAttribute(QgsField("acq_datetime_PST", QVariant.String))

        layer.updateFields()

        # Populate the new field

        for feature in layer.getFeatures():

            feature["acq_datetime_PST"] = feature["acq_datetime_PST"]

            layer.updateFeature(feature)

        layer.commitChanges()

        print_message("Applying custom styling to the layer...")

        ranges = []

        for class_def in [

            {"name": "< 6h", "min": 0, "max": 6, "color": "purple"},

            {"name": "6 - 12h", "min": 6, "max": 12, "color": "red"},

            {"name": "12 - 24h", "min": 12, "max": 24, "color": "orange"},

            {"name": "24 - 48h", "min": 24, "max": 48, "color": "yellow"},

            {"name": "48h - 7days", "min": 48, "max": float("inf"), "color": QColor(204, 255, 0)}

        ]:

            symbol = QgsSymbol.defaultSymbol(layer.geometryType())

            symbol.setColor(QColor(class_def["color"]))

            symbol.setSize(12)

            symbol.symbolLayer(0).setStrokeWidth(1.2)

            range = QgsRendererRange(class_def["min"], class_def["max"], symbol, class_def["name"])

            ranges.append(range)

        renderer = QgsGraduatedSymbolRenderer('time_difference_hours', ranges)

        renderer.setMode(QgsGraduatedSymbolRenderer.Custom)

        layer.setRenderer(renderer)

        layer.triggerRepaint()

        print_message("Preparing KML exports...")

        kml_output_dir = os.path.expanduser(f"~/desktop/FIRMShotspots_scriptRun_{timestamp}")

        os.makedirs(kml_output_dir, exist_ok=True)

        class_definitions = [

            {"name": "less_than_6h", "min": 0, "max": 6},

            {"name": "6_to_12h", "min": 6, "max": 12},

            {"name": "12_to_24h", "min": 12, "max": 24},

            {"name": "24_to_48h", "min": 24, "max": 48},

            {"name": "48h_to_7days", "min": 48, "max": float("inf")}

        ]

        for class_def in class_definitions:

            print_message(f"Exporting KML for hotspots {class_def['name']}...")

            class_layer = QgsVectorLayer("Point?crs=EPSG:4326", f"Fires {class_def['name']}", "memory")

            class_provider = class_layer.dataProvider()

            class_provider.addAttributes(layer.fields())

            class_layer.updateFields()

            features = []

            for feature in layer.getFeatures():

                time_diff = feature['time_difference_hours']

                if class_def['min'] <= time_diff < class_def['max']:

                    features.append(feature)

            class_provider.addFeatures(features)

            class_kml_path = os.path.join(kml_output_dir, f"fires_{class_def['name']}.kml")

            save_options = QgsVectorFileWriter.SaveVectorOptions()

            save_options.driverName = "KML"

            QgsVectorFileWriter.writeAsVectorFormatV2(class_layer, class_kml_path, QgsCoordinateTransformContext(), save_options)

        print_message("Exporting map extent to KML...")

        extent_layer = QgsVectorLayer("Polygon?crs=EPSG:4326", "Clipping_Extent", "memory")

        extent_provider = extent_layer.dataProvider()

        extent_provider.addAttributes([QgsField("id", QVariant.Int)])

        extent_layer.updateFields()

        extent_feature = QgsFeature()

        extent_feature.setGeometry(QgsGeometry.fromRect(extent))

        extent_feature.setAttributes([1])

        extent_provider.addFeature(extent_feature)

        extent_kml_file_path = os.path.join(kml_output_dir, "clipping_extent.kml")

        QgsVectorFileWriter.writeAsVectorFormat(extent_layer, extent_kml_file_path, "utf-8", extent_layer.crs(), "KML")

        print_message("\n" + "*" * 80)

        print_message("*" * 80)

        print_message(f"*     KML FILES HAVE BEEN SAVED TO:     *")

        print_message(f"*     {kml_output_dir}     *")

        print_message("*" * 80)

        print_message("*" * 80 + "\n")

print_message("Script completed successfully. FIRMS hotspot data has been processed and loaded into QGIS.")