How do I create blank geotiff with same spatial properties as existing geotiff?

Question

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.)

Answer 1

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

Answer 2

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

Answer 3

The GDAL version 3.3 (in year 2021) introduced a new option for the utility gdal_create.

-if <input_dataset> New in version 3.3.

Name of GDAL input dataset that serves as a template for default values of options -outsize, -bands, -ot, -a_srs, -a_ullr and -a_nodata. Note that the pixel values will not be copied.

Documentation includes an usage example

Initialize a blank GeoTIFF file from an input one:

gdal_create -if prototype.tif output.tif

Answer 4

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).

Answer 5

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

Answer 6

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!