Generate high-resolution tiled imagery from USDA NAIP data
I use OpenMapTiles to create self-hosted vector map tiles. However, that doesn't provide aerial imagery.
The US Department of Agriculture (USDA) captures high-resolution aerial imagery for the continental US. The USGS's web portal allows for easy downloading of NAIP imagery.
According to the NAIP website:
NAIP imagery is acquired at a one-meter ground sample distance (GSD) with a horizontal accuracy that matches within six meters of photo-identifiable ground control points, which are used during image inspection.
The style JSON spec tells Mapbox GL how to style your map. Add the raster imagery tiles as a source to overlay them with the other sources.
Within sources, each object key defines the name by which the later parts of
style.json should refer to the layer. Note the difference between a normal
raster layer and the terrain RGB layer.
"sources": {
"openmaptiles": {
"type": "vector",
"url": "https://api.maptiler.com/tiles/v3/tiles.json?key={key}"
},
"naip": {
"type": "raster",
"url": "https://example.com/url/to/tile.json",
"tileSize": 512
}
}Where the tile.json for a raster layer should be something like:
{
"attribution": "<a href=\"https://www.usgs.gov/\" target=\"_blank\">© USGS</a>",
"description": "NAIP Imagery",
"format": "png",
"id": "naip",
"maxzoom": 15,
"minzoom": 11,
"name": "naip",
"scheme": "tms",
"tiles": ["https://example.com/url/to/tiles/{z}/{x}/{y}.png"],
"version": "2.2.0"
}Later in the style JSON, refer to the raster to style it. This example shows the raster layer between zooms 11 and 15 (inclusive), and sets the opacity to 0.2 at zoom 11 and 1 at zoom 15, with a gradual ramp in between.
{
"id": "naip",
"type": "raster",
"source": "naip",
"minzoom": 11,
"maxzoom": 15,
"paint": {
"raster-opacity": {
"base": 1.5,
"stops": [
[
11,
0.2
],
[
15,
1
]
]
}
}
}Clone the repository:
git clone https://github.com/nst-guide/naip
cd naip
This is written to work with Python >= 3.6. To install dependencies:
pip install -r requirements.txt
This also has dependencies on some C/C++ libraries. If you have issues installing with pip, try Conda:
conda env create -f environment.yml
source activate naip
Downloads USGS elevation data for a given bounding box.
> python download.py --help
Usage: download.py [OPTIONS] [FILE]...
Download raw NAIP imagery for given geometry
Options:
--bbox TEXT Bounding box to download data for. Should be
west, south, east, north.
-b, --buffer-dist FLOAT Buffer to use around provided geometry. Only
used with --file argument.
--buffer-unit [mile|meter|kilometer]
Units for buffer. [default: mile]
--buffer-projection INTEGER EPSG code for projection used when creating
buffer. Coordinates must be in meters.
[default: 3488]
--overwrite Re-download and overwrite existing files.
--help Show this message and exit.
This script calls the National Map API and finds all the 3.75'x3.75' NAIP imagery files that intersect the given bounding box or geometry. By default, it only downloads the most recent image, if more than one exist.
The script then downloads each of these files to data/raw/. By default,
it doesn't re-download and overwrite a file that already exists. If you wish to
overwrite an existing file, use --overwrite.
Use gdalbuildvrt to generate a virtual dataset of all image tiles and
gdal2tiles to cut the output raster into map tiles.
gdal2tiles.py options:
--processes: number of individual processes to use for generating the base tiles. Change this to a suitable number for your computer.- I also use my forked copy of
gdal2tiles.pyin order to generate high-res retina tiles
First, download desired DEM tiles, unzip them, build a VRT (Virtual Dataset),
and optionally download my fork of gdal2tiles which allows for creating
512x512 pngs.
# Download for some geometry
python download.py -b 2 example.geojson
# Create virtual raster:
# The -b flags take each of the RGB bands at full opacity.
# The -addalpha sets areas with no source raster as transparent.
# See https://github.com/nst-guide/naip/issues/1
gdalbuildvrt -b 1 -b 2 -b 3 -addalpha data/naip.vrt data/raw/*.jp2
# Download my fork of gdal2tiles.py
# I use my own gdal2tiles.py fork for retina 2x 512x512 tiles
git clone https://github.com/nst-guide/gdal2tiles
cp gdal2tiles/gdal2tiles.py ./
# Generate tiled imagery
# --exclude excludes transparent tiles from output tileset
# You may want to set max zoom level 15
./gdal2tiles.py --processes 10 --exclude data/naip.vrt data/naip_tilesRunning this for a .5 mile buffer for the entire Pacific Crest Trail, with 512x512 pixel tiles and a max zoom of 16 generated ~60GB of png files. Around ~45GB of these tiles are in just zoom level 16, and while I'll probably set a max zoom of 15 in the future, ~15GB is still more than I want to serve.
Note that these CLI commands worked on macOS, and some might need different syntax on Linux.
Note: It appears from limited inspection that Mapbox GL JS shows WebP images slightly pixelated? See #2.
Google's webp format is a great compressor of png images. Running the cwebp
cli (you can download from Google's site or through homebrew ) and setting the
quality to 80/100 reduces file size by about 85%
> cwebp 10299.png -q 80 -o 10299.webp
By default this is lossy compression. You can set -lossless if you'd like,
but lossy is fine for my needs.
The docs for the -q option say:
Specify the compression factor for RGB channels between 0 and 100. The default is 75.
In case of lossy compression (default), a small factor produces a smaller file with lower quality. Best quality is achieved by using a value of 100.
In case of lossless compression (specified by the -lossless option), a small factor enables faster compression speed, but produces a larger file. Maximum compression is achieved by using a value of 100.
I think this is really confusing because a higher -q value is less compressed
when running lossy and more compressed when running lossless!
To create a full hierarchy of webp images:
cd data/naip_tiles
for f in */*/*.png; do mkdir -p ../naip_tiles_webp/$(dirname $f); cwebp $f -q 80 -o ../naip_tiles_webp/$f ; doneUnfortunately, webp isn't supported everywhere, namely on older browsers and on all iOS browsers.
Because WebP images aren't supported on all devices, I also need to serve a PNG layer. Unlike the terrain-rgb elevation layer, for which lossy compression would erase all meaning of the encoded elevation values, for regular images lossy compression won't make much of a visual difference.
From limited searching, the best open source, command line PNG compressor seems
to be pngquant. This creates PNG images that can be
displayed on any device. You can install with brew install pngquant.
In order to create a directory of compressed images:
cd data/naip_tiles
find . -name '*.png' -print0 | xargs -0 -P8 -L1 pngquant --ext -comp.png --quality=70-80Setting quality to 70-80 appears to create files that are about 25% of the
original size.
I couldn't figure out an easy way to write the output png files to a new directory. By default they're written in the same directory as the original file, with an added extension.
# in data/naip_tiles already
mkdir ../naip_tiles_png
# First move all tiles with extension -comp.png to the new directory
# Take out --remove-source-files if you want to copy, not move the files
# I use rsync because I couldn't figure out an easy way with `mv` to move while
# keeping the directory structure.
rsync -a --remove-source-files --include "*/" --include="*-comp.png" --exclude="*" . ../naip_tiles_png/
# Then rename all of the files, taking off the -comp suffix
cd ../naip_tiles_png
find . -type f -name "*-comp.png" -exec rename -s '-comp.png' '.png' {} +When run on Halfmile's PCT track data with a .5 mile buffer (including alternates), the resulting data after pngquant compression has file sizes:
> du -csh * | sort -h
| Directory size | Zoom level |
|---|---|
| 8.0K | 4 |
| 20K | 5 |
| 36K | 6 |
| 108K | 7 |
| 348K | 8 |
| 1.2M | 9 |
| 4.1M | 10 |
| 15M | 11 |
| 60M | 12 |
| 240M | 13 |
| 979M | 14 |
| 3.8G | 15 |
| 15G | 16 |
| 20G | total |
I'll probably serve up to zoom 15 through the web browser, and maybe up to zoom 14 for offline download.