# gdalwarp DEM with different resolution extents

Geographic Information Systems Asked by TazMainiac on December 5, 2020

I’m trying to convert DEM source (LIDAR is original source, in ASC format) to a custom internal format that has specific resolution and alignment values. I also want to generate different output resolutions (each lower resolution is 1/9th resolution of higher res).

I’ve written a GDAL plugin to convert from source DEM to the internal format, but it requires the source data to already be in the correct resolution and alignment. So I’m trying to figure out the right combination of gdalwarp, gdal_translate and gdal_retile to achieve this. I need the pixel spacing (resolution) to match the target format resolution and I need those pixels to align with the target formats restrictions. I’m fine cutting data off the edges of the source.

I use python to calculate the correct commands and run them. I end up with (target SRS is WGS84/EGM96 – hopefully I got the t_srs right).

gdalwarp -t_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=./egm96_15.gtx" -tap -et 0 -r bilinear -tr <res> <res> -te <res_xmin> <res_ymin> <res_xmax> <res_ymax> <src.vrt> <dest.tif>


where <res> is specific to the resolution (such as 0.00001028806584362140 for 1/27th AS). The -te values are specific to the alignment requirements for that specific resolution. So I figure out the coordinates of the target format file that is "inside" the source data (cutting off the edge), and make sure the extent is a multiple of the target file size (again staying inside the bounds of the source data).

Then I use gdal_retile to split up <dest.tif> into individual files to fit what my custom format needs (300×300, 1 pixel overlap, so I do 301×301).

Finally I use gdal_translate with my custom plugin to convert the split files to my custom format.

This seems to work for my different output resolutions, but when I try and view the data, the higher res data does not line up with the lower res data – when I load the higher res data over the top of the higher res data, terrain features do not line up as expected.

So first question: when gdalwarp process -te does it try and fit the source data to that extent (i.e. ‘squashing’ the source data to fit, assuming -te is smaller than source extent), or does it ‘cut off’ the part of the source data that is outside that extent? I’m trying to get it to cut it off. If it’s ‘squashing’ instead, that might be why I’m seeing alignment problems.
[edit : comment indicates that gdalwarp "cuts" the data]

I’ve also tried:

gdalwarp -t_srs "..." -tr <res> <res> <src.vrt> <temp.tif>
gdal_translate <res_xoff> <res_yoff> <res_xsize> <res_ysize> <temp.tif> <dest.tif>


i.e. just have gdalwap reproject without trying to "trim" or "align" the output, then use gdal_translate to do the trimming to match my desired alignment. But this seems to end up with the same problem, and I’m concerned that the pixel alignment won’t match up.

Does this sound like I’m going about this the right way? Any suggestions to try and figure how to re-project and trim the data to match resolution and alignment restrictions?

Thanks!

gdalwarp will reproject your data while preserving the geographic properties of your input data. In other words, it will be as close as possible from the initial information at each geographic location (this information will be slightly altered if the position of the output pixel is not exactly the same as the position of the pixel in the output file).

-te defines the size of the output file in the coordinate system of the output file (by default, the coordinate system of the output file = the coordinate system of the input). If larger than the input, the extent will be filled with NoData values. If smaller than the input, the input will be cut and some of the information is "lost" You need to provide the extent as geographic information.

-tr allows you to define the output size horizontally and vertically. If you change the grid, gdal will use the infomation of the input pixels that are close to the center of the new pixel to compute a new value (this is called resampling, managed with the -r parameter, by default the value of the nearest pixel is used)

get information from gdalinfo to determine the exact value because of potential rounding issues, no need to use -tap if you have the precise -te and -tr. Use the infomation from gdal and not from other software, because the pixel coordinate used for the extent could vary dependng of the software (e.g. center or top left of the pixel)

for the second step (changing resolution), in your case the best way is to use gdal_translate because, once you have the correct extent, you can control the number of pixels and this will operate a resampling

for example, if your input file is 27000 * 27000 pixel, you can divide the number of pixels by 9 by reqesting an output file of 3000 * 3000 pixels without changinf the extent.

gdal_translate -outsize 3000 3000 temp.tif dest.tif

As a remark, 301 cannot be divided by 9, so your pixels will not overlap correctly.

## Related Questions

### QGIS modelbuilder, having a choicebox/value map of strings as string input

1  Asked on November 18, 2020

0  Asked on November 18, 2020 by bryan

### Using mouse-over on raster surface to stream values to text file using ArcGIS Desktop?

1  Asked on November 18, 2020 by gisi

### Iteratively polygonize MultiLineString

0  Asked on November 18, 2020

### How to assess the visibility of a feature from anywhere on the landscape (DEM)

0  Asked on November 17, 2020 by argosy

### Joining attributes by location for more than one located feature using R

2  Asked on November 17, 2020

### Pattern did not match any band error for Sentinel 2, 2015 year image in Google Earth Engine

1  Asked on November 17, 2020 by arushi

### Edit GISJOIN attribute ID from census data

0  Asked on November 16, 2020 by user65127

### How can I export images with their metadata from Google Earth Engine to Google Cloud Storage?

1  Asked on November 16, 2020 by saurabh-khanna

### Creating point precisely on line

2  Asked on November 16, 2020 by sweet-sugar-cola

### Trying to export shapefile attributetable to JSON using GeoPandas without the geometry

2  Asked on November 16, 2020 by lord-corwin

### Assign ICEWS events to administrative regions from GDAM data

0  Asked on November 15, 2020 by tl_bayern

### Lock error when converting polygon to raster in ArcGIS Desktop?

2  Asked on November 15, 2020 by tony-bof

### Digitized contours on map and interpolate contour line to create surface

1  Asked on November 14, 2020 by fanfangis

### Load and apply Style from Postgres in QGIS 3

2  Asked on November 13, 2020 by simu

### Two PopupTemplates

1  Asked on November 12, 2020 by pfalbaum

### GeoPackage Lookup Tables

0  Asked on November 12, 2020 by narmaps