Counting trees from the air, and proving the coordinates are right

I built tkmedia-gis-tree-detection to do two things: detect tree crowns on a georeferenced aerial tile, then write each detection out as a real-world coordinate. On one NEON test tile it finds 55 crowns and exports them as a GeoJSON point layer. The repo is MIT licensed.

The detection step is standard. Coordinate handling decides whether the output is useful, because a mistake there still produces a file that opens and draws in QGIS, only in the wrong place.

What the pipeline does

The pipeline reads a georeferenced GeoTIFF, detects crowns, converts each box centre to a real-world coordinate and writes a GIS-ready point layer.

The code follows this sequence:

1. rasterio.open() on the input tile, pulling crs, transform, and pixel size off the dataset. Two asserts guard the two quiet failures that can be caught at read time: a missing CRS, and an identity transform. A GeoTIFF with an identity transform is a plain image wearing a GeoTIFF extension, and every coordinate you derive from it will be wrong in a way that looks fine.
2. DeepForest with the pretrained weecology/deepforest-tree weights, pinned to CPU.
3. predict_image() when the raster fits one patch, predict_tile() when it does not. The sample tile is 400×400 px, so it takes the first path.
4. Filter boxes below the score threshold, default 0.3.
5. Box centre to coordinate, then reproject.
6. Write GeoJSON in EPSG:4326, GeoPackage in the native CRS, an annotated PNG, and a summary.json.

The sample tile is NEON airborne imagery over the Ordway-Swisher Biological Station in Florida, at 0.1 m ground sample distance, native CRS EPSG:32617. Detection takes a few seconds on a laptop CPU.

The transform, and the offset argument inside it

The box centre goes through the raster's own affine transform. Three lines do the important part:

col_c = (xmin + xmax) / 2.0
row_c = (ymin + ymax) / 2.0
easting, northing = transform_xy(transform, rows=row_c, cols=col_c, offset="ul")

transform_xy is rasterio.transform.xy. The affine comes from the file, not from a value I supplied. The georeferencing already lives in the GeoTIFF; my job is to carry it through unchanged.

offset="ul" tells rasterio to treat the row/col you hand it as measured from the upper-left corner of the pixel grid rather than from pixel centres. Pass offset="center" instead and every point shifts by half a pixel on both axes. At 0.1 m GSD that is 5 cm, which will never look wrong on a map and will quietly disagree with any other layer you overlay.

Then to_crs(epsg=4326) for the GeoJSON, because GeoJSON is defined in WGS 84; UTM values written into a .geojson produce points thousands of kilometres from the site.

How the coordinates were verified

I checked them by hand against the affine coefficients read from the file.

results/VERIFICATION.md recomputes one detection:

easting  = c + a * col + b * row = 404211.9 + 0.1 * 351.93216 = 404247.0932
northing = f + d * col + e * row = 3285142.9 + (-0.1) * 367.17075 = 3285106.1829

The value the script wrote for that tree is 404247.093, 3285106.183. Agreement to three decimal places, which is where the rounding stops. The document also records the offset="center" alternative and the 0.05 m per-axis shift it produces, so the choice is written down rather than assumed, and a tile-centre sanity check that lands at −81.989891, 29.692504, in Florida, where OSBS is.

There is no automated test. No tests/, no CI, no assertion in the code that the transform is right. The verification is a markdown file a human wrote once. The arithmetic is reproducible, but nothing will notice the day someone changes offset and reruns. It is the first thing I would fix before anyone depended on this pipeline.

What comes out

55 points, each carrying score, easting, northing, native_crs, and a grid_zone of NW, NE, SW or SE. Per-zone counts are 10, 15, 14 and 16.

I split the count by quadrant because "how many trees" is rarely the real question. "How many in this parcel" is. The zones are a crude stand-in for a real parcel boundary, but they make the output useful: a count you can group by area instead of one number you have to trust.

The highest score in the set is 0.7987.

Where it breaks

The 55 is a detection count at threshold 0.3. It is not ground truth. Nobody walked that plot and counted trees, so the correct reading is "the model found 55 crowns it was at least 30% confident about", not "there are 55 trees".

Beyond that:

• The weights are a generalist pretrained detector with no fine-tuning on this site. Move to a different canopy type, a different season, or a different sensor and the score distribution moves with it.
• Positional accuracy is the crown-box centroid pushed through the raster's georeferencing. It inherits whatever error the orthophoto already has.
• It reports the crown centre, not the trunk. A leaning tree on a slope puts the trunk somewhere the box centre never sees. For forestry inventory that gap matters and this pipeline does not model it.
• Single tile, CPU, one patch. The tiled path exists but the sample never exercises it.

The repo, including the verification document, is at github.com/SM1LE21/tkmedia-gis-tree-detection.

I build this kind of pipeline on a client's own orthophotos and parcel boundaries, with the coordinate checks above written in as tests. If you have imagery to count or map, write to me.

Imagery: NEON Airborne Observation Platform, Ordway-Swisher Biological Station, distributed with the open-source DeepForest sample data.