top of page
logo
Code Archive

Australia Is on Fire, According to a Satellite with a Clipboard — Part I

  • Writer: Mic
    Mic
  • Jun 23
  • 10 min read

Updated: 2 days ago

A satellite fire dataset looks wonderfully authoritative until you remember what it actually contains: points. This is not an exact fire perimeter or incident boundaries. You should not expect a neat polygon saying “everything inside here burned”.

NASA FIRMS — the Fire Information for Resource Management System — gives us active fire and thermal anomaly detections at satellite overpass time. That is extremely useful, but it is also the kind of useful that needs a warning label.

In plain Python terms, FIRMS lets us ask NASA for rows where a satellite saw something hot enough to look fire-like. The rows have longitude, latitude, acquisition time, confidence, and a few useful measurement fields. From there we can build maps, timelines, and spatial summaries.

In geospatial terms, however, the first lesson is more cautious: a FIRMS row is an observation, not an incident report. The satellite gives us a hot pixel. It does not give us the full fire history, the official perimeter, the local warning zone, the cause, the fuel conditions, or the emotional state of the nearest eucalyptus tree.

The danger here is that: the data looks map-ready, but it needs interpretation before we start telling stories with it.

This first part has one job: build the clean, cached dataset and create a first overview map of the selected fire case locations across Australia. We will not yet create any graphs or plots, but focus on acquiring and reading the data. The only thing we will create visually is a plot that shows where everything takes place.

In the second part we then use the cleaned data to make individual fire maps, daily detection curves, Fire Radiative Power timelines, and hotspot-density grids.

I guess ChatGPT thinks that satellites are cool...

The tools doing the unglamorous work

The updated script uses the usual Python geospatial suspects, plus contextily for the basemaps used in Part II.

So our main packages will be

  • Pandas

  • Requests

  • GeoPandas

  • Shapely

  • Matplotlib

  • Contextily

The roles are straightforward. We will use requests to talk to the FIRMS API, then pandas to read and clean the CSV data, use geopandas to turn longitude and latitude columns into geometry, use shapely to create points and bounding boxes, and finally matplotlib to draw the figures. In some of them we will use contextily to add web tiles underneath the projected maps later in the workflow.

In other words, nobody is doing anything mystical. The satellite does the orbital drama. Python mostly does tidy columns, polite requests, and rectangles with delusions of geographic precision.

The script starts with imports, local folders, API constants, and plotting defaults.

There are a few choices hidden in that setup.

MAX_DAY_RANGE = 5 is there because the FIRMS Area API is queried in short date windows. Instead of asking for a large historical period all at once, the script walks through each case study in blocks of up to five days. That makes the requests predictable and keeps the code easier to debug when one chunk fails.

DATA_DIR = Path("data/firms") gives the downloaded CSV files a proper home. This matters because maps are usually adjusted many times. You do not want to re-download the same NASA data every time you move a title, resize a marker, or decide that one shade of red is not quite apocalyptic enough.

OUTPUT_DIR = Path("outputs") is where the final images go. The script creates it immediately so the plotting functions can save figures without each one worrying about folder housekeeping.

The source is set to VIIRS_SNPP_SP, which means the VIIRS instrument on Suomi-NPP using the standard processing product. For historical examples, that is a sensible default because standard processing products are meant for stable retrospective use. For near-real-time monitoring, check FIRMS data availability first and choose one of the near-real-time products instead.

BASEMAP_SOURCE = cx.providers.CartoDB.VoyagerNoLabels matters in the next part. The detailed fire maps use a background tile map with as few labels as possible. Labels are useful when navigating, but they can also fight with fire markers, legends, and annotation boxes. In the attached image examples we used cx.providers.OpenStreetMap.Mapnik instead to make the colours stand out a bit more, but for a more analytic map, the labels are a bit too much clutter. In those cases, the basemap is there to provide geographic context, not to shout town names over the data.


Choosing fire cases without worshipping the rectangle

We will be looking at eight Australian fire case from the last few years. Each case is stored as a small dataclass.

We have the full name for titles, slug for filenames, and label for...well map labels. The short label will be mostly useful for the plot at the end of this post.

The case list looks like this:

Two of these are from a summer in Australia, that is usually called the Black Summer. We will see in Part II why that name is quite warranted in these two cases. Each case has six fields, the first three we already explained above. Afterwards, bbox is the FIRMS Area API bounding box in the order 

west, south, east, north.

This order is worth repeating because bounding boxes are where perfectly nice scripts go to quietly become wrong. It is not north/south/east/west, and it is not latitude first. Longitude west, latitude south, longitude east, latitude north. Then finally start and end define the inclusive query window.

The boxes are intentionally rectangular and blunt. A box around Kangaroo Island includes water. A box around a mainland case may include nearby activity that is not part of the named incident. That is not a bug in the code. It is the trade-off of using a rectangular API query for real-world fire behaviour, which tends to have poor respect for neat rectangles.

We have seen in the GIS series in part I and II how to cut this into geometric regions if ones need a more specific observation, like for a specific council area.


Asking FIRMS in small date chunks

FIRMS does not need us to dramatically ask for everything at once. The script is better behaved if it asks in small chunks, caches the answers, and then joins them together locally.

The Area API URL has a regular structure: map key, source product, bounding box, day range, and start date. We will create the url in the following function.

The line doing most of the work is this one:

area = ",".join(f"{value:.4f}" for value in bbox)

It transforms a tuple like this:

(146.0, -38.4, 150.4, -36.0)

into the coordinate string expected by the API:

146.0000,-38.4000,150.4000,-36.0000

Then the script creates date chunks:

The generator starts at the case-study start date, moves forward in blocks of up to five days, and yields a pair like ("2019-12-20", 5). FIRMS interprets that as the start date plus the requested inclusive day range.

The final chunk of a case may be shorter than five days, so the function yields the actual day_range rather than blindly returning 5 every time. Without that, the script would ask for dates beyond the intended case window. The satellite may be orbiting Earth, but our query should not drift past the end date because we got lazy with loop boundaries.


Turning API replies into a local dataset

The next problem pretty boring: external APIs do not always fail in the tidy way we would choose for a tutorial. Sometimes the response is empty. Sometimes it says there were no fires. Sometimes it is an error message with the emotional tone of a locked filing cabinet.

So the script reads each FIRMS CSV response defensively.

There are three useful safety checks here.

First, timeout=60 stops the request from hanging. Then, response.raise_for_status() turns HTTP problems into visible errors. Finally, the column check makes sure the response is actually a fire-detection table before the rest of the workflow starts treating it like one.

Now the downloading, caching, and chunking come together in one function.

If the cache file already exists and refresh=False, the function reads the CSV from disk instead of calling NASA again. If the cache is missing, it loops over the date chunks, downloads each piece, collects non-empty DataFrames, combines them with pd.concat, drops duplicates, cleans the result, and writes the cleaned file to disk.

The cache filename includes both the case slug and the source product:

cache_path = DATA_DIR / f"{case.slug}_{source}.csv"

That means a future version can compare VIIRS and MODIS for the same case without overwriting its own files. The tiny sleep(0.2) is not really interesting, but it is polite. This is not be overlooked. We are talking to a public service, keep requests to a minimum amount needed and give it ample time.


The timestamp goblin

FIRMS acquisition time arrives as a compact value such as 57, 1432, or 2359. That is convenient for CSV storage and less convenient for humans who want a timestamp.

The cleaning function fixes the problem and also adds a few case-study columns.

The quiet hero is this line:

acq_time = df["acq_time"].astype(str).str.zfill(4)

It turns 57 into 0057, which can be read as 00:57. After that, pd.to_datetime combines the date and time into a proper UTC timestamp.

The function also creates acq_day, because daily summaries are much easier when the date part has already been extracted. Finally, the script maps VIIRS confidence codes to readable labels and keeps only nominal and high confidence detections:

df = df[df["confidence"].isin(["n", "h"])].reset_index(drop=True)

This is a pragmatic filter. It does not make the dataset perfect, but it avoids building the plots around low-confidence detections. The result is still not a fire perimeter, but at least we are not asking every suspicious warm pixel to audition for a disaster movie.


From rows to geometry

So far, the data is a table. A useful table, but still a table. To make maps, GeoPandas needs a geometry column.

The important method is:

gpd.points_from_xy(df["longitude"], df["latitude"])

Notice the order, here points are built from x first, y second, which means longitude first and latitude second. This is another geospatial detail that rewards obsessive repetition. In everyday usage, latitude usually comes first, but most geometry libraries expect x/y.

The CRS is set to EPSG:4326. This means, the usual WGS84 longitude/latitude coordinate system. Later, when we add web basemaps, the data will be projected to Web Mercator. For the Australia overview map in this post, longitude and latitude are enough.


Making figures that explain themselves

The individual plots in this project are meant to live inside a blog post. Without much of a report around them, we will add a small interpretation guide. A title and legend are good. A small explanation box is better.

The key detail is transform=ax.transAxes. It means the text position is expressed relative to the axis instead of in longitude/latitude or Web Mercator coordinates. A y-position like -0.20 places the annotation just below the plot panel.

The function also adjusts the bottom margin so the box is not cut off when the image is saved. This is not geospatially deep, but it prevents a very common plotting annoyance: the figure looks fine in the notebook, then the exported PNG quietly amputates the caption.

Static plots should carry their own legend, title, and short interpretation guide wherever possible. Otherwise, they become beautifully exported riddles.


A first map needs case centres, not every fire point

Before plotting thousands of detections, it helps to answer a simpler question: where are the selected cases?

For that, the overview map does not need every FIRMS point. It needs one marker per case study, placed at the centre of that case’s bounding box.

These centre points are not ignition points, incident centroids, or official locations. They are simply centre points of the rectangular query areas. That distinction is important enough to put directly into the figure. The marker says “the query box is around here”, not “the fire politely started here”.


A first map of the case locations

Now we are ready to map out where the locations are. This is the first genuinely map-like output of the project, but it is intentionally modest. It gives the reader a mental scaffold before Part II starts showing individual fire maps, timelines, and density grids.

The country outline comes from Natural Earth:

world = gpd.read_file(NATURAL_EARTH_COUNTRIES)
australia = world[world["NAME"] == "Australia"].copy()

Then the function plots Australia in a quiet background style and overlays the case centre points in crimson. The labels use the short label field from FireCase, so the map stays readable even with eight cases.

The label background colour is intentionally light red:

"facecolor": "mistyrose"

That makes the labels feel visually connected to the fire markers without turning the whole map into a warning sign. The alpha value keeps the boxes semi-transparent, and the edge is removed so they do not become heavy little postage stamps.

The explanation box below the map prevents the most obvious misreading:

Red markers = centre points of the bounding boxes
used for each fire case
These are overview locations only

A centre marker is just a centre marker. It is not the fire front, it is not the worst affected point, and it is definitely not the satellite’s secret favourite pixel.


Running the first half of the workflow

To run the first part, fetch the cases, combine the non-empty frames, print the detection counts, and call the overview map function.

After this runs, the cache folder should contain one CSV per case and source product, and the output folder should contain:

outputs/firms_fire_case_locations_australia.png

The printed counts are useful as a sanity check. They tell you whether data actually arrived for each case and whether one query window is much larger than the others. They are not, by themselves, a scientific comparison of fire severity. Detection counts depend on area, duration, satellite overpass timing, cloud and smoke, fuels, fire intensity, and the bounding box chosen by the person currently pretending their rectangle is a research instrument...that would be me. The resulting plot should look something like this.

First plot of this project, simply marking what we look at.

The workflow so far

At this point, the script has built the unglamorous part of the project:

FireCase definitions
        ↓
FIRMS Area API URLs
        ↓
Five-day date chunks
        ↓
CSV responses
        ↓
Cleaned and cached DataFrames
        ↓
GeoPandas point helpers
        ↓
Australia-wide case overview map

That is a good place to stop Part I because the next stage changes the question. Here we asked: can we get the data into a reliable local form, and where are the selected cases?

Part II asks a different set of questions. Where did detections cluster inside each case box? When did activity spike? How does Fire Radiative Power change over time? What happens when thousands of points are aggregated into a simple grid?

In short: Part I turns the satellite clipboard into a tidy table and a first map. Part II starts waving the coloured markers around.

The complete source for this part as well as the next can be found in this GitHub repository.


Sources and further reading

In contrast to other posts, we are using a few outside API's and data. Hence some sources for further reading:

Comments

bottom of page