3

Using GDAL, how can I simply create an empty single band raster? I can't seem to find out how. I know in PostGIS you can create a new table with a 'raster' column type, then insert a band into it. However I'm looking for a slick, GDAL solution (or another efficient non-PostGIS solution).

What I'm looking to create is a raster (GeoTIFF) with these settings:

format: -of GTiff

type: -ot UInt16 (16 bit unsigned integer)

cell resolution: -tr 5 5 (5x5m)

nodata: -a_nodata 255 (the values I intend to burn into this raster range from 1 to 10)

SRID: 27700

extent: -te 0 0 680000 1240000 (rectangular extent of Great Britain)

resulting raster would contain almost 34 billion cells (at 5x5m cell size)!

Improve this question
7
  • 1
    You're asking two questions - how can I create an empty raster?, then how can I populate it?. Can you edit your question to contain just a single question. In this case I would suggest not trying to create an empty raster, but just generating a distance or proximity raster directly from the vector data. Though you may need to tile it and mosaic it together at the end, that's a big raster.
    – user2856
    Commented Jan 24, 2020 at 12:18
  • thanks. Edited. Regarding size, sure the coverage is large but file size on a small sample area of 13690x13690m (2738x2378 pixels) is under 1mb. Extrapolating that out to the whole country, it should be a manageable size of a couple of Gbs.
    – Theo F
    Commented Jan 24, 2020 at 12:21
  • You are setting 255 to no data value, while this would be unhandy if you have 2^16 numbers available to you. Perhaps you mean to set 2^16-1 as the no data value?
    – Johan
    Commented Jan 24, 2020 at 13:02
  • @Johan hi. My values will eventually only range from 1 to 10, so 255 will never be reached anyway.
    – Theo F
    Commented Jan 24, 2020 at 13:33
  • 1
    Yes, byte and unsigned 8 bit hold 0-255
    – Johan
    Commented Jan 24, 2020 at 13:50

2 Answers 2

Reset to default
4

In Python, this would be:

import gdal
import osr

driver = gdal.GetDriverByName('GTiff')

spatref = osr.SpatialReference()
spatref.ImportFromEPSG(27700)
wkt = spatref.ExportToWkt()

outfn = '/path/to/out.tif'
nbands = 1
nodata = 255
xres = 5
yres = -5

xmin = 0
xmax = 680000 
ymin = 0
ymax = 1240000 
dtype = gdal.GDT_Int16

xsize = abs(int((xmax - xmin) / xres))
ysize = abs(int((ymax - ymin) / yres))

ds = driver.Create(outfn, xsize, ysize, nbands, dtype, options=['COMPRESS=LZW', 'TILED=YES'])
ds.SetProjection(wkt)
ds.SetGeoTransform([xmin, xres, 0, ymax, 0, yres])
ds.GetRasterBand(1).Fill(0)
ds.GetRasterBand(1).SetNoDataValue(nodata)
ds.FlushCache()
ds = None
Improve this answer
6
  • thanks however that appears flipped. ie. it appears above my area of interest.
    – Theo F
    Commented Jan 24, 2020 at 14:11
  • on ds.SetGeoTransform replace ymax with ymin to prevent the raster being vertically flipped into the map region above your intended area.
    – Theo F
    Commented Jan 24, 2020 at 14:49
  • 1
    Thanks for the sharp eye. From the documentation on the gdal raster data model: The (GT(0),GT(3)) position is the top left corner of the top left pixel of the raster. So, the ymax should be given as I put in my example. What should be done to flip the image to what you want, the yres should be negative. I updated the example so that it works now.
    – Johan
    Commented Jan 24, 2020 at 15:21
  • 1
    In general, If you imaging a plot, the top left corner coordinate of the plot will be (xmin, ymax). So, if you want the image to go down (normally southward), you should "add" negative pixels, i.e. subtract y resolution. To check myself, I looked up some geotransforms of some tiffs that I could get my hands on: they all have a negative y-resolution. So i think that the documentation that you saw has it right.
    – Johan
    Commented Jan 24, 2020 at 16:20
  • 1
    @TheoF The documentations are correct. The "origin" for most image processing applications is the top-left rather than bottom-left. Johan is correct that this is accounted for by having a negative y resolution, which does seem odd.
    – Jon
    Commented Jan 24, 2020 at 16:21
4

Since version 3.2, gdal offers gdal_create command for this purpose. Just simply use the following:

gdal_create -of GTiff -ot UInt16 -a_nodata 255 -burn 0 -outsize 680000 1240000 -a_srs "EPSG:27700" -co COMPRESS=LZW large_file.tif

I assumed by the extent, you meant the width and height of the image. For the cell resolution of 5m or whatever else, just set up your georeferenced bounds of the output file based on the image size and the pixel resolution. Add this option to the command line: -a_ullr <ulx> <uly> <lrx> <lry>

I also added the compression option of -co COMPRESS=LZW to be able to save disk space for this huge empty file! However, it did not work! It seems to be a bug for such a big image size. The compression works for smaller image size.

Also note the -burn value has been set to 0. You can change it to whatever value you like. It takes a while for the image to be created; my write speed was 0.4 MB/s so it took around 2 hours and 6 minutes and the file size was ~2.5 GB. Be patient!

Source: https://gdal.org/programs/gdal_create.html

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.