5

How do I use gdal to create a blank geotiff with the same spatial properties (size, resolution, extent, projection) as an existing geotiff?

(The reason that I want to do this is that I'm using gdal_rasterize to rasterize a polygon shapefile to overlay with existing geotiffs, and I think it will be simpler to burn the values into a matching layer rather than specifying all the properties in the rasterize command.)

  • Are you wanting the extent to be different to the extent of your vector data? – MappaGnosis Dec 8 '16 at 14:14
  • It can be the same extent. – user1521655 Dec 8 '16 at 14:24
  • In which case, this process is an unnecessary step and you don't need to specify a heap of extent and SRS properties - see my answer below. – MappaGnosis Dec 8 '16 at 15:49
4

One line solution: modify the last gdal_calc.py example on http://www.gdal.org/gdal_calc.html:

gdal_calc -A input.tif --outfile=empty.tif --calc "A*0" --NoDataValue=0

Check the result:

gdalinfo empty.tif -hist
Driver: GTiff/GeoTIFF
Files: empty.tif
Size is 10, 10
Coordinate System is `'
Origin = (950.000000000000000,1050.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (     950.000,    1050.000)
Lower Left  (     950.000,      50.000)
Upper Right (    1950.000,    1050.000)
Lower Right (    1950.000,      50.000)
Center      (    1450.000,     550.000)
Band 1 Block=10x10 Type=Byte, ColorInterp=Gray
0...10...20...30...40...50...60...70...80...90...100 - done.
  256 buckets from -0.5 to 255.5:
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  NoData Value=0
  • For large rasters this is very slow (even using such options like --co SPARSE_OK=TRUE), presumably because it needs to read through the entire A raster (not just the metadata). Isn't there any other way that will work almost instantaneously? – benjimin Jan 22 '18 at 5:57
  • @benjimin, gdal_merge.py with the -createonly option seems to work instantaneously. I will write another answer soon, please have a look at it. – user30184 Jan 22 '18 at 11:05
  • This alternative (using gdal_translate to rescale to zero) is much faster than gdal_calc but not instantaneous: gdal_translate -ot Byte -of GTiff -scale 0 255 0 0 -co TILED=YES -co COMPRESS=LZW -co BIGTIFF=YES -co SPARSE_OK=TRUE input.tif empty.tif – benjimin Jan 22 '18 at 22:39
  • It is still writing a value for each pixel. I know that GDAL can somehow initialize an image without filling the pixels. The result is not useful for other software than GDAL but for gdal_rasterize it should do. – user30184 Jan 23 '18 at 5:52
3

One more way to do it with GDAL is to use gdal_merge.py http://www.gdal.org/gdal_merge.html with -createonly and -init options.

Usage:

gdal_merge -createonly -init "0 0 0" -o empty.tif sentinel.tif
1

Instead, I think that you need to specify all the properties in the rasterize command in order to obtain a correct result. As example, if need to rasterize the sample_shape.shp shapefile and you want to use the sample_raster.tif raster as a reference for the new output, you may use its properties to create your output (called output_raster):

from osgeo import gdal, ogr

# Filename of input OGR file
sample_vector = 'sample_shape.shp'

# Filename of the raster Tiff that will be created
sample_raster = 'sample_raster.tif'

# Filename of the raster Tiff that will be created
output_raster = 'output_raster.tif'

# Define pixel_size and NoData value of new raster
pixelSizeX = sample_raster.rasterUnitsPerPixelX()
pixelSizeY = sample_raster.rasterUnitsPerPixelY()
NoData_value = -9999

# Open the data source and read in the extent
source_ds = ogr.Open(sample_vector)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# Create the destination data source
x_res = int((x_max - x_min) / pixelSizeX)
y_res = int((y_max - y_min) / pixelSizeY)
target_ds = gdal.GetDriverByName('GTiff').Create(output_raster, x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixelSizeX, 0, y_max, 0, -pixelSizeY))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

# Rasterize
gdal.RasterizeLayer(output_raster, [1], source_layer, burn_values=[0])

(This code is adapted from Convert an OGR File to a Raster).

  • What language is this? – user1521655 Dec 8 '16 at 14:28
  • PyQGIS, i.e. Python scripting for QGIS – mgri Dec 8 '16 at 14:31
  • Any idea how to do this with the command-line utility for Linux/bash? – user1521655 Dec 8 '16 at 14:33
  • No, I'm sorry, but I think this could be another question... However, if you use QGIS, it should be easy to adapt this code to your needs. – mgri Dec 8 '16 at 14:36
0

As you can use the same extent as your vector file, this process is unnecessary. You don't need to create the empty file. As of GDAL 1.8 and above, it will take the extent of the vector file. Check out the documentation for te. Likewise the project will match the vector file too. You can use the -init option to initialize your background raster value to avoid nodata too.

You only need to specify the resolution but that is less work than creating the empty raster in the first place!

  • But I would still need to specify resolution, right? – user1521655 Dec 8 '16 at 15:54
  • Specifying the resolution is less typing than it takes to create the empty raster. – MappaGnosis Dec 8 '16 at 15:59
  • Not if I want to automate this with scripts and can't count on always wanting the same output resolution. – user1521655 Dec 8 '16 at 16:00
  • You just use a couple of variables in the script for x and y. It is very straightforward and you'd still have the same problem with creating your empty raster first too. – MappaGnosis Dec 8 '16 at 16:01
0

One option is to copy the image with gdal_translate and set the scale to use value range from zero to zero.

Usage example

gdal_translate -scale 0 400000 0 0 sentinel.tif empty.tif

Gdalinfo with histogram (gdalinfo image.tif -hist) from source and target. Remember to delete the old histogram data first by removing the *.aux.xml file if that exists.

Source:

  Metadata:
    STATISTICS_MAXIMUM=18119
    STATISTICS_MEAN=3572.9272485136
    STATISTICS_MINIMUM=187
    STATISTICS_STDDEV=1787.6104593137

Target:

  Metadata:
    STATISTICS_MAXIMUM=0
    STATISTICS_MEAN=0
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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