How do I create blank geotiff with same spatial properties as existing geotiff?
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
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
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).