Cutting large DEMs into square, partially overlapping tiles

I’ve been doing a lot of shaded relief using Blender, and one practice I’ve developed is feeding Blender square DEM tiles.

The DEMs I use in maps tend to be too big for my installation of Blender to handle. One I had recently was 6348 x 10899 pixels; another was 4155 x 3572.  So that’s one reason for tiling them: I can produce shaded relief of the tiles separately in Blender and then mosaic the results.

And as long as I’m tiling, I’d like to use square tiles. That’s because when you subdivide a plane in Blender you have to use the same factor on both X and Y axes. Call me superstitious, but I get uncomfortable knowing that I’ve mapped 1803 pixels of DEM onto 1500 cells along one axis, and 2403 pixels to the same 1500 cells on the other. It seems to me this should cause loss of precision.

By feeding Blender square tiles, I eliminate that. I can subdivide the plane 2048 times on both axes, and feed it a 2050 x 2050 DEM.

This is not the kind of tiling that you would use a tool like gdal2tiles.py to do. That tool produces tiles that do not overlap, but I want some overlap. This is for a couple of reasons.

One is that Blender sometimes produces a bit of an edge artifact on a tile’s down-light edge (that is, the edge of the tile away from the light source). I can trim this off, and if there’s a bit of overlap between adjacent tiles, I get a nice smooth mosaic.

Left: mosaic of two non-overlapping tiles -- note the "seam" down the centre. Right: the left tile trimmed, and the right tile (with an overlap of 200 pixels) laid under it. Light is from 337.5°

Left: mosaic of two non-overlapping tiles — note the “seam” down the centre. Right: the left tile trimmed, and the right tile (with an overlap of 200 pixels) laid under it. Light is from 337.5°

Another reason for overlap is that if a shadow-casting feature (a mountain, a cliff) is just off the up-light edge of a tile, Blender can’t know to cast the shadow from that feature onto the tile. Overlap cleans this up by putting a healthy margin of features on the up-light side.

How much you overlap seems to depend on the topography and the scale. With small scale mapping, you don’t need much overlap because the mountains can’t cast their shadows very far. At large scales they can cast them many, many pixels.

Comparison of shaded relief produced for 1:1,000,000 (left), 1:150,000 (centre), and 1:30,000 (right).

Comparison of shaded relief produced for 1:1,000,000 (left), 1:150,000 (centre), and 1:30,000 (right).

To make tiling-with-overlap easier, I’ve written a bash script, called maketiles.sh, to tile large images, which you can download. This script uses the gdal_translate tool, so you also need to have GDAL installed.

Usage looks like this (assuming linux here):

./maketiles.sh <image file> <desired tile size in pixels> <number of columns in image> <number of rows in image> <tile_filename_prefix> <overlap>

It’s the rare DEM that can be cut perfectly into tiles without the last row and column being a bit off-sized. This script resolves that by increasing overlap on the final row and column so that all tiles are the requested size.

Let’s run through an example.

I have a DEM called 093L03-4-5-6_dem_26909.tif which is 4155 x 3572 pixels.

(If I don’t initially know its dimensions, I can get them through gdalinfo:

gdalinfo 093L03-4-5-6_dem_26909.tif
Driver: GTiff/GeoTIFF
Files: 093L03-4-5-6_dem_26909.tif
093L03-4-5-6_dem_26909.tfw
Size is 4155, 3572)

Let’s say I want tiles that are 2000×2000 pixels with a 100 pixel overlap, and I want the resulting tiles (and their world files) to be named tile_1_1.tif, tile_1_2.tif, etc.

The command would be

./maketiles.sh 093L03-4-5-6_dem_26909.tif 2000  4155 3572 tile 100

The script first echos back your parameters, and then tells you how many tiles will be produced. In this case we’ll be getting 6 tiles in 2 rows and 3 columns. It also tells you how much extra overlap you’ll have on the last row and column.

image file is 093L03-4-5-6_dem_26909.tif
desired tile size is  2000 x 2000
image is  4155 columns by  3572  rows
prefix is tile
overlap is 100 pixels.

Estimating 3 columns of tiles (X), by 2 rows of tiles (Y).

Taking a closer look at columns.
When laying out 3 columns of tiles, it looks like the final column will fall at 5800 pixels.
That's bigger than the image, so we'll back off the final column of tiles to begin at 2155 and end exactly at 4155.
This final column of tiles will overlap the second-to-last column by 1745 pixels.

Taking a closer look at rows.
When laying out 2 row of tiles, it looks like the final row will fall at 3900 pixels.
That's bigger than the image, so we'll back off the final row of tiles to begin at 1572 and end exactly at 3572.
This final row of tiles will overlap the second-to-last row by 428 pixels.

Layout has 6 tiles in 2 rows and 3 columns.

Finally, it computes the maximum overlap you could have without increasing the number of tiles, and asks if you want to try that.

Would you like to go with an overlap of 428, which would spread the overlap out but keep the same number of tiles? [Y/n]

If you accept this increased overlap, it confirms the change, and then you get a final chance to approve the process or exit.

Overlap set to 428.
overlap is 428.

Hit Enter to proceed with cutting tiles, or Ctrl-C to exit.

Assuming you hit Enter/Return, it then calls gdal_translate many times as is necessary, giving you extent info for each tile.

Producing tile tile_1_1.tif (row 1, column 1).  Extents: (0, 0) to (2000, 2000).
gdal_translate  -co "TFW=YES" -srcwin 0 0 2000 2000 093L03-4-5-6_dem_26909.tif tile_1_1.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.

Producing tile tile_1_2.tif (row 1, column 2).  Extents: (1572, 0) to (3572, 2000).
gdal_translate  -co "TFW=YES" -srcwin 1572 0 2000 2000 093L03-4-5-6_dem_26909.tif tile_1_2.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.

The final tile in the row will be a bit too big.
To keep all tiles the same size, I'm backing off the upper left corner of this tile from (3144,0) to (2155,0).

Producing tile tile_1_3.tif (row 1, column 3).  Extents: (2155, 0) to (4155, 2000).
gdal_translate  -co "TFW=YES" -srcwin 2155 0 2000 2000 093L03-4-5-6_dem_26909.tif tile_1_3.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.
ROW IS FINISHED.

Producing tile tile_2_1.tif (row 2, column 1).  Extents: (0, 1572) to (2000, 3572).
gdal_translate  -co "TFW=YES" -srcwin 0 1572 2000 2000 093L03-4-5-6_dem_26909.tif tile_2_1.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.

Producing tile tile_2_2.tif (row 2, column 2).  Extents: (1572, 1572) to (3572, 3572).
gdal_translate  -co "TFW=YES" -srcwin 1572 1572 2000 2000 093L03-4-5-6_dem_26909.tif tile_2_2.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.

The final tile in the row will be a bit too big.
To keep all tiles the same size, I'm backing off the upper left corner of this tile from (3144,1572) to (2155,1572).

Producing tile tile_2_3.tif (row 2, column 3).  Extents: (2155, 1572) to (4155, 3572).
gdal_translate  -co "TFW=YES" -srcwin 2155 1572 2000 2000 093L03-4-5-6_dem_26909.tif tile_2_3.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.
ROW IS FINISHED.
IMAGE IS FINISHED.

Tiles are named with the prefix you gave, followed by _X_Y.tif, where X is the row position, and Y the column position, of the tile. So in this case we got tile_1_1.tif, tile_1_2.tif, etc., along with tile_1_1.tfw, tile_1_2.tfw, etc

Advertisements

5 thoughts on “Cutting large DEMs into square, partially overlapping tiles

  1. Pingback: Shaded relief with BlenderGIS, part 1 | The Wandering Cartographer

  2. Hello,

    First of all, thanks a lot for sharing your script here! I think it can be very useful. I was wondering whether you could help me with a technical problem… I’m running your script in Windows in MSys. It runs beautifully, however at the end it displays ROW IS FINISHED IMAGE IS FINISHED. However there are no newely created files, I would expect them in the original folder where the input tif is.
    Do you have any idea where it could have gone awry?

    Greetings from the Netherlands,

    Jeroen

    Like

    • Hi Jeroen,

      Wow, I had no idea one could run a linux script under Windows!

      If you read the script, you’ll see that the line that creates the output files is a call to the GDAL utility “gdal_translate”, such as

      gdal_translate -co “TFW=YES” -srcwin $presentXMin $presentYMin $tileSize $tileSize $imageFile $tileFileName

      Do you have GDAL installed? Can you test it from the command line, using your image file and some suitable parameters?

      Like

      • Hello Morgan,

        Thanks for the reply! Yes GDAL is installed, it was installed throught the OSGeo4W utility, and I believe that it is being called correctly in your application as it goes through all the steps required. It is just the output that is not being produced. Any ideas?

        Like

      • Hi Jeroen,

        Two ideas:

        1) Is the image that you are tiling a geoTiff? If the TIF image does not have embedded projection information in it (or an accompanying TFW file), gdal_translate won’t know what to do with it.

        2) Is your geoTiff projected in a coordinate system whose units are metres? (E.g., UTM or Gauss-Kruger). If, for example, you have an image that is in decimal degrees, you will first have to convert it to a local metric projection.

        If both of these are true, I would test gdal_translate outside the script. A line like

        gdal_translate -co “TFW=YES” -srcwin 0 0 100 100 myImage.tif outputImage.tif

        should give you a 100×100 image result. It will wind up in the directory from which you run the script.

        – Morgan

        Like

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