top of page
Search
Writer's pictureMarco G J

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

Updated: Oct 22

(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.")

14 views0 comments

Comments


bottom of page