geotiff.js: How to Get Projected GeoTIFF Data for a Latitude-Longitude Coordinate
Reprojecting latitude and longitude coordinates into a GeoTIFF’s coordinate system using Javascript, and querying data with geotiff.js
Javascript has been used for many years to provide interactive webmaps that are usually composed of vector data and RGB image tiles. These simplistic frontends have historically required more capable languages and server technologies to actually provide and render the geospatial data being visualized, in addition to mechanisms for querying or analyze them.
As the Javascript language and ecosystem has matured, what was once almost impossible has become straightforward and even decently performant. Thanks to the work of dedicated open source developers, it is now fairly easy to work directly with geospatial raster data both in the frontend (browser) or backend (NodeJS). There is even multithreading support!
Let’s take a look at one piece of functionality that used to be quite difficult to do in pure Javascript: querying geospatial raster data for a specific coordinate. These datasets are most typically available in GeoTIFF format, and to work with that we’ll use one of the more popular libraries for working with the file format: geotiff.js.
It is easy to use this library to extract pixel data from a GeoTIFF. But, given a latitude-longitude coordinate, how do you figure out what pixel to read from? What if the GeoTIFF is in a projected coordinate system?
Unlike GDAL (a popular opensource geospatial library, which has bindings for NodeJS), geotiff.js has no reprojection functionality: it merely parses GeoTIFF files. This means our solution is going to be a little more involved than the equivalent approach in GDAL, which, at its simplest, is a single command:
gdallocationinfo -valonly -wgs84 <geotiff filename> <longitude> <latitude>We’ll just need to write some code to mirror the extra work that GDAL is doing under the hood: extracting the projection information from the raster, reprojecting the latitude-longitude coordinates into the raster’s coordinate system, and determining which pixel of the raster that point applies to.
Why Not Just Reproject the Entire GeoTIFF First?
This entire process can mostly be circumvented by reprojecting our GeoTIFF into a latitude-longitude coordinate system like EPSG:4326. But for many projected coordinate systems, this introduces serious inaccuracies.
For example, let’s consider querying the NBM weather model to get the forecasted temperature for the summit of Little Bear Peak, Colorado:
This example was chosen as it is near the edge of a pixel in the dataset.
Now let’s reproject the dataset into EPSG:4326 using nearest-neighbor resampling. The same coordinate now returns a different value:
We can improve the reprojection slightly by using a different resampling method, such as bilinear:
But clearly, we’ve lost some of the original weather model information. We can reduce how much information is lost by making the reprojected raster much, much larger, but even at massive sizes there will still be edge cases where locations near the boundary of a pixel will not pull the original weather model data correctly.
Whether or not this is a big deal depends on the purpose and resolution of your data. If you’re using a high resolution digital elevation model to construct elevation profile graphs for routes that are miles long, then it’s probably not much of a concern if the elevation data is slightly off in places. Conversely, if you’re using weather model data that has coarse resolution, then being in the wrong grid cell could result in a significantly different forecast for an area.
Overall: it’s best to use the original projection for querying raw data from a geospatial raster.
Reprojecting Latitude and Longitude to the Raster Coordinate System
The first step is to get the GeoTIFF’s projection information and use that to reproject our desired coordinate(s) into the raster’s coordinate system.
A large part of what makes a GeoTIFF a geospatial raster, instead of a regular TIFF image, is the use of GeoKeys (per the OGC GeoTIFF spec). We can read the GeoKeys for a given raster using geotiff.js:
// assume the variable gtiff is the geotiff object
// that was created by loading data with geotiff.js
const image = await gtiff.getImage();
const geoKeys = image.getGeoKeys();We will pass these GeoKeys to a projection library which will handle the heavy lifting, as projections tend to be quite complicated!
For instance, weather model data for North America tends to be in a Lambert conformal conic projection. If we start with a latitude-longitude coordinate derived from an ellipsoidal datum (commonly WGS84), then we’re faced with a rather complicated formula:
No need to reinvent the wheel, as the open source community has poured years of works into comprehensive reprojection libraries which easily handle all sorts of transformations between complicated map projections.
For Javascript, one of the most robust, popular, and full-featured libraries is Proj4js, which is an offshoot of the common proj library that powers GDAL. We’ll just need to generate a proj string using our raster’s GeoKeys and pass the final string to the library.
However, there are a lot of valid GeoKeys:
/**
* Geokeys. If you're working with `geotiff` library, this is result of `image.getGeoKeys()`.
* @typedef {Object} module:geokeysToProj4.GeoKeys
* @property {number} GeographicTypeGeoKey See GeoTIFF docs for more information
* @property {number} GeogGeodeticDatumGeoKey See GeoTIFF docs for more information
* @property {number} GeogPrimeMeridianGeoKey See GeoTIFF docs for more information
* @property {number} GeogLinearUnitsGeoKey See GeoTIFF docs for more information
* @property {number} GeogLinearUnitSizeGeoKey See GeoTIFF docs for more information
* @property {number} GeogAngularUnitsGeoKey See GeoTIFF docs for more information
* @property {number} GeogAngularUnitSizeGeoKey See GeoTIFF docs for more information
* @property {number} GeogEllipsoidGeoKey See GeoTIFF docs for more information
* @property {number} GeogSemiMajorAxisGeoKey See GeoTIFF docs for more information
* @property {number} GeogSemiMinorAxisGeoKey See GeoTIFF docs for more information
* @property {number} GeogInvFlatteningGeoKey See GeoTIFF docs for more information
* @property {number} GeogPrimeMeridianLongGeoKey See GeoTIFF docs for more information
* @property {number} ProjectedCSTypeGeoKey See GeoTIFF docs for more information
* @property {number} ProjectionGeoKey See GeoTIFF docs for more information
* @property {number} ProjCoordTransGeoKey See GeoTIFF docs for more information
* @property {number} ProjLinearUnitsGeoKey See GeoTIFF docs for more information
* @property {number} ProjLinearUnitSizeGeoKey See GeoTIFF docs for more information
* @property {number} ProjStdParallel1GeoKey See GeoTIFF docs for more information
* @property {number} ProjStdParallel2GeoKey See GeoTIFF docs for more information
* @property {number} ProjNatOriginLongGeoKey See GeoTIFF docs for more information
* @property {number} ProjNatOriginLatGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseEastingGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseNorthingGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseOriginLongGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseOriginLatGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseOriginEastingGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseOriginNorthingGeoKey See GeoTIFF docs for more information
* @property {number} ProjCenterLongGeoKey See GeoTIFF docs for more information
* @property {number} ProjCenterLatGeoKey See GeoTIFF docs for more information
* @property {number} ProjCenterEastingGeoKey See GeoTIFF docs for more information
* @property {number} ProjCenterNorthingGeoKey See GeoTIFF docs for more information
* @property {number} ProjScaleAtNatOriginGeoKey See GeoTIFF docs for more information
* @property {number} ProjScaleAtCenterGeoKey See GeoTIFF docs for more information
* @property {number} ProjAzimuthAngleGeoKey See GeoTIFF docs for more information
* @property {number} ProjStraightVertPoleLongGeoKey See GeoTIFF docs for more information
* @property {number[]} GeogTOWGS84GeoKey Datum to WGS transformation parameters, unofficial key
*/It may be worth using a helper library like geotiff-geokeys-to-proj4 to handle the proper translation of geotiff.js’s geokeys object into a proj string. In fact, the above comment block is from that library.
With those two libraries installed, we don’t need much to code to set up our reprojection functionality:
import proj4 from 'proj4';
import * as geokeysToProj4 from 'geotiff-geokeys-to-proj4';
... // Not shown: importing geotiff.js and loading a geotiff file.
// for this example, the variable `gtiff` references the file.
const image = await gtiff.getImage();
const geoKeys = image.getGeoKeys();
const projObj = geokeysToProj4.toProj4( geoKeys );
const projection = proj4( `WGS84`, projObj.proj4 );We’ve now created a proj projection object called “projection.” It’s very easy to use:
const { x, y } = projection.forward( {
x: -105, // the longitude
y: 40 // the latitude
} );The above code reprojects the latitude-longitude coordinates into the raster coordinate system (likely meters or feet).
Now we need to get the actual raster pixel associated with those new coordinates. First, let’s get some information about the size of the raster and its bounding box:
const width = image.getWidth();
const height = image.getHeight();
const [ originX, originY ] = image.getOrigin();
const [ xSize, ySize ] = image.getResolution();
const uWidth = xSize * width;
const uHeight = ySize * height;The above code establishes a relationship between the size of the raster in pixels and the actual coordinates/size of the raster in the projected coordinate system.
We then take our projected coordinate relative to the origin coordinates of the raster and see where it falls in relation to the width and height of the raster.
// x and y come from the projection.forward example earlier
const percentX = ( x - originX ) / uWidth;
const percentY = ( y - originY ) / uHeight;
const pixelX = Math.floor( width * percentX );
const pixelY = Math.floor( height * percentY );Now we have a pixel value for the original latitude-longitude coordinate!
You’ll probably want to ensure this pixel is actually inside the raster. If the original coordinate didn’t intersect the raster, then pixelX or pixelY could be negative or greater than the width/height of the raster.
Extracting the GeoTIFF Data for the Pixel
If the pixel is inside the raster, then we have two methods of extracting the data, considering a GeoTIFF with one single band.
The first method is to pass a 1x1 pixel window to geotiff.js’s “readRasters” function:
const [ value ] = await image.readRasters( {
interleave: true,
window: [ pixelX, pixelY, pixelX + 1, pixelY + 1],
samples: [ 0 ]
} );This method is straightforward and obvious, but there is a small overhead to calling “readRasters” which could add up if you are querying thousands of coordinates.
In that case, there may be a slight performance benefit with a different approach — at the cost of RAM, which could be considerable for very large GeoTIFFs. We can read the entire raster into memory as an ArrayBuffer:
const data = await image.readRasters( {
interleave: true,
samples: [ 0 ]
} );Then we just need to grab the corresponding index:
const value = data[ width * pixelY + pixelX ];In practice, geotiff.js caches tiles and stripes internally, so the performance gain is only slim at best, and generally unnoticeable for modest datasets.
Thus, this approach will only provide a (minor) speed benefit if querying a lot of data that is spread across a majority of the stripes or tiles in the GeoTIFF. Otherwise, the overhead of reading the entire image, including tiles or stripes that are unused, will outweigh the potential performance gain, though this can be mitigated by reading a window that is sized to the minimum and maximum x and y or the coordinate collection.
Using the above techniques, it is possible to implement the following features in an interactive web map, either entirely in the browser or with a basic NodeJS API:
- Clicking on a web map to get raster data for a specific point
- Labelling point features with underlying raster data (e.g. a weather map with city names showing the forecast temperature)
- Performing basic geometry operations against raster data (such as an intersection operation using a collection of points or a bounding box), and displaying that information in a table or exporting it as a CSV
- Generating charts or lists of data for a 2d cross-section or line (for instance, an elevation profile)
- Performing statistical analysis and data transformations on raw GeoTIFF data
This sort of functionality (let alone actually rendering geospatial raster data onto a map) has traditionally required specialized GIS server software and dedicated infrastructure, with interactive web maps doing little more than sending requests to the server and displaying the response. As the Javascript ecosystem evolves, browsers are becoming more capable of parsing geospatial data and thus reducing reliance on specialized backend technology.
Traditionally, a popular method to make GeoTIFF/geospatial raster data available to web maps was having a GIS server application providing common Open Geospatial Consortium interfaces — such as WMS (Web Map Service) or WMTS (Web Map Tile Service). These interfaces would allow the browser to request image files or tiles of the data, which would be rendered server-side and sent to the browser as RGB image files.
The server could also provide a WCS (Web Coverage Service) interface, which would let the browser query the underlying geospatial raster data for specific coordinates. If using proprietary GIS server technology (such as Esri’s ArcGIS products), the endpoints and syntax may be a little different, but the end result is the same.
However, the recent maturity of the cloud-optimized GeoTIFF standard plus the ability to access and render the data directly in the browser using openlayers + geotiff.js has removed the need for a GIS server providing WMS and WMTS endpoints, at least for some scenarios. The querying techniques described in this article show one method of removing the need for a WCS or query endpoint.
Instead, one can merely host COGs using their existing cloud provider (for instance, an S3 bucket) and make the web browser responsible for parsing and working with the geospatial data. Organizations with NodeJS backends can consolidate GIS functionality into these stacks and reduce the need for specialized software and additional hardware.
Finally, another huge benefit of being able to do more in the browser with geospatial raster data is the reduction or elimination of the need for internet connectivity. For smartphone mapping applications, particularly those intended for use in remote areas, this provides an avenue to provide offline geospatial tools for apps that are already written in Javascript.
Several of these apps are very popular and provide functionality such as slope angle calculation, terrain profile charts, saved weather maps, etc. These apps will let you save maps or areas for offline usage — but once the user is offline, many features stop working as they rely on servers to actually parse and analyze the geospatial data.
These features don’t have to break once internet access is lost — the technique mentioned in this article is one of many that can allow for the app to remain functional while offline.
