Interpolation and resampling across projections and spatial resolutions with GDAL

It is very common in the physical sciences to work with data streams that have different spatial resolutions and sometimes even with different native projections. Hence, the need emerges to resample or interpolate data from one resolution to another. In the following example I take PRISM climatological precipitation data, nominally at 1/24 degree (2.5 arc-minutes, ~4 km) horizontal resolution and resample it to a 150 m equal area UTM projection.

The target projection is for a river basin in southwestern Colorado called the Uncompahgre R, with an area of approximately 386 km2.

The raw PRISM data (Available here: http://www.prism.oregonstate.edu/) are available in ASCII ArcInfo, or ASCII grid format, i.e. with a 6 row header, over the conterminous U.S. (CONUS) domain.

  1. The first step is optional, namely clipping the raw file to a smaller area that is easier to work with. This step is uses the ‘gdal_translate’ utility with the “-projwin” option, specifying the upper-left and lower-right corners of the bounding box.

gdal_translate -projwin -115 45 -100 35 <InFile.asc> <ClippedFile.asc>

Notice the above coordinates represent the maximum (45) and minimum (35) latitude, and the maximum (-100) and minimum (-115) longitude.

  1. Next, we’ll convert the clipped data to “.tif” format, which will allow for processing with the gdal tools. “.tif” is the default format when converting formats with gdal_translate. It should be noted that:

gdal_translate <ClippedFile.asc> <ClippedFile.tif>

  1. Now the key step of remapping to the 150 m UTM grid with gdalwarp, using a bilinear method, being sure to assign the correct geo-referencing:

gdalwarp -s_srs WGS84 -t_srs ‘+proj=utm +zone=13 +datum=WGS84’ <ClippedFile.tf> <ClippedFileResampled.tif> -r bilinear -tr 150 150 -tap

A few key options are that we assign the source data (latitude longitude grid) a ‘WGS84′ projection, and the destination a UTM’ projections, for zone 13. The resampling method is blinear and the target grid size is 150m by 150m. Note a full list of UTM zones can be found here: http://www.dmap.co.uk/utmworld.htm

  1. If you wish to clip the resampled file to a smaller target window (in this case the boundary of the Uncompahgre R), then you’ll need to first convert a grid for that target domain to .tif format (using gdal_translate above) and then use the ‘gdalinfo’ command to get the extents, e.g:

gdal_translate <TargetMask.asc> <TargetMask.tif>

gdalinfo <TargetMask.tif>

Now using the upper-left and lower-right corners of the target mask, clip the resampled file to the target grid:

gdalwarp -te 238236.485 4180122.372 291336.485 4247322.372 <ClippedFileResampled.tif> <TargetResampledFile.tif>

I included those numbers for max/min coordinates above just to show that these are actually in units of meters(!) for UTM projections.

  1. Lastly, you may wish to convert these data either back to ASCII, or else to netCDF, e.g.

ASCII:

gdal_translate -of AAIGrid -co DECIMAL_PRECISION=2 <TargetResampledFile.tif> <TargetResampledFile.asc>

Notice, the DECIMAL_PRECISION of the output data can be specified.

netCDF

gdal_translate -of netCDF <TargetResampledFile.tif> <TargetResampledFile.nc>

Advertisements

3 thoughts on “Interpolation and resampling across projections and spatial resolutions with GDAL

  1. Thanks for posting Ben, this is really helpful. Quick question — do you know of an easy way to extract the data value at a particular lat/lon? Either from the command line or in a language-specific wrapper.

    For Windows people, the GDAL libraries come bundled with the Python(x,y) installation (https://code.google.com/p/pythonxy/) which will let you do things like `gdalwarp` from the Cygwin terminal. (Or the Windows cmd but, eh).

    Of course on a more developer-friendly OS you can `brew install gdal` or `sudo apt-get install gdal-bin`.

  2. Hi Jon,
    I am not immediately aware of how to do this, other than once the data are in netCDF (step 6), then you can use ncks -H -v -d,lat,[value or range] -d,lon,[value or range] [infile.nc]

  3. Hey I want to interpolate point data into a raster file using gdal python …..code u provide me with code and some explanations because I m a cs student and have very little knowledge of all this !!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s