Nav Icon

Back Arrow Postholer.Com
Resource for Hikers

Log in or register for full access to site features & free downloads

Printed Maps & App
Stay found. Stay informead. Our premium printed topographic maps and matching app for a once in a lifetime adventure

OnLine Data Books
Free online data books for your favorite trail, distances, elevation, climate, way points, terrain & fauna

Interactive Maps
GPS aware interactive maps, your location, topo, hybrid, satellite, trail track, way points, road, weather/snow overlays

Active Fires & Smoke
Know before you go! Get the lastest active fire & smoke information for your favorite trail

Snow Conditions
Get the lastest snow conditions for your favorite trail. Quantity, coverage, SNODAS, MODIS, historical and multi-year comparisons

Trip Planners
Create an extensive hike plan easily configurable to your hiking style. Distances, days, resupply, access points, etc

Gear Builder
Create your own fully customizable gear list with weight, pricing and divisable by section.

Trail Animals
Exhaustive resource for species ranges along your favorite trail. Includes amphibians, birds, mammals and reptiles.

Complete Gear Lists
Hundreds of complete gear lists by those that walked the walk.

Elevations Profiles
The big picture. Single elevation profile of your favorite trail.

Journal Tools
Search for a journal, create a journal, add/edit an entry, configure your journal, EMail updates, integrated interactve trail map, PLB locations and more

Wall Maps
Print out your favorite trail to 6 feet high. Elevation chart and resupply locations.

Postholer Forum
Source for trail and site information or just talk about your favorite trail

Where it all begins.

Faster Raster Data Acquisition

July 11, 2023

femafhz map
Double click on map for production example

To the point...
How to retrieve pixel values from 15 different rasters at a given coordinate in 1 line of code.

...and a bit more verbose
At some point you'll want your interactive map users to retrieve raster value(s) at a given map click. This is often done with a web service on a back end server/service. How you do this can be extremely complex or insanely simple. It all depends on the tools you use.

Will you do this with an ArcPy script or AGOL or maybe a python script and all the attendant environment? Think about this for a moment. What does your environment require to retrieve a single pixel value from a raster using ESRI or Python? What resources does the care and feeding of that environment require? I've been there and the answer makes me cringe.

The short answer to my rhetorical question.
We've used gdalbuildvrt, specifying paths to all 15 rasters, to create a virtual raster, femafhz.vrt. We'll use Bash and GDAL to get 15 pixel values from 15 different rasters at a single coordinate:


   gdallocationinfo -geoloc -valonly /opt/gis/vrt/femafhz.vrt -100.123456 40.987654 | \
   awk -v div=15 'NR%div { printf("%s,", $0); next} {print $0}'
The above will return CSV with 15 different values:
Great. Wonderful. What about performance? What about performance at scale? Well, let's time it:
   time ./
   0.098u 0.009s 0:00.10 90.0%	0+0k 0+0io 0pf+0w
It took 1/10 of a second to return 15 raster pixel values. That's reasonable. BUT, what if we have many coordinates. Surely, we'll see scale issues. Let's find out. The VRT file I'm using covers CONUS. I've generated 100,000 random CONUS coordinates sorted by longitude and placed them in 'coords.txt'. We'll modify our script to accomodate these values:

   gdallocationinfo -geoloc -valonly /opt/gis/vrt/femafhz.vrt < coords.txt | \
   awk -v div=15 'NR%div { printf("%s,", $0); next} {print $0}'
With 100K random coordinates it took 31.78 seconds or 3,146 coords per second. Put another way, at 15 rasters per coordinate, that's 47,199 raster pixel values per second. Is that adquate? You tell me.
   time ./
   28.953u 2.583s 0:31.78 99.2%	0+0k 0+12936io 0pf+0w
These tests were performed on a production AWS EC2 small instance with 2GB ram and 2.5Ghz processor. On a quiet server you'll get even better times. Our VRT file contains 15 rasters with pixel resolutions between 10 and 30 meters. The VRT file was created using gdalbuildvrt command.

Part of the irony here is ESRI uses GDAL under the hood for format conversions. Any geo-spatial in Python will require you to install GDAL libraries. Why not just use GDAL directly where you can and cut out the middleman?

Postholer.Com © 2005-2024 - Sitemap - W3C - @postholer - GIS Portfolio