Create and save raster dataset using GDAL in Python

GDAL can not only read, but also create data sets. There are two ways for GDAL to create a dataset: one with the Create() method and the other with the CreateCopy() method. Which method should be used depends on the data and on the other hand, depending on the format of the file. All drivers that support the creation of new files support the CreateCopy() method, but only the Create() method is supported in part. You can check the value of DCAP_CREATE and DCAP_CREATECOPY in the meta information of the driver to determine whether the Create() or CreateCopy() method is supported. As shown in the following code:

>>> from osgeo import gdal
>>> format = "GTiff"
>>> driver = gdal.GetDriverByName( format )
>>> metadata = driver.GetMetadata()
>>> if gdal.DCAP_CREATE in metadata and metadata[gdal.DCAP_CREATE] == 'YES':
...     print ('Driver %s supports Create() method.' % format)
...
Driver GTiff supports Create() method.
>>> if gdal.DCAP_CREATECOPY in metadata and metadata[gdal.DCAP_CREATECOPY] == 'YES':
...     print ('Driver %s supports CreateCopy() method.' % format)
...
Driver GTiff supports CreateCopy() method.

In the above example, first create a GeoTiff dataset format driver driver, and then extract the driver metadata to see if there is a DCAP_CREATE or DCAP_CREATECOPY attribute in the metadata. By judging the values of these two properties, you know that the driver supports Create() and CreateCopy() methods.

Not every data supports the Create() method, such as JPEG, PNG and other formats are not supported. There is a trick for these formats to write out the data format in memory. You can use Create() to create an intermediate file in GeoTiff format, and then use the CreateCopy() method to create JPEG or PNG.

Create images using the CreateCopy function

The CreateCopy() function is relatively simple to use, which is to convert an image of one format directly into an image of another format, which is equivalent to a format conversion.

Here’s an example of creating an image using the CreateCopy() function:

>>> import gdal
>>> src_filename = "/gdata/geotiff_file.tif"
>>> dst_filename = "/tmp/xx_geotiff_copy.tif"
>>> driver = gdal.GetDriverByName('GTiff')
>>> src_ds = gdal.Open( src_filename )
>>> dst_ds = driver.CreateCopy( dst_filename, src_ds, 0 )

The first argument to the CreateCopy() function is the source file name, and the second function is the target file name. The third and subsequent parameters are optional. If not entered, the program runs by default. The value of the third parameter is 0 or 1, and the value is 0. Even if the original data cannot be accurately converted to the target data, the program is still executed. CreateCopy()method, no fatal error will occur. This kind of error may be that the output format does not support the data type of the input data format pixel, or the target data does not support the write spatial reference and so on. This function also has a fourth parameter, and the fourth parameter refers to some special parameters that are used in the format conversion.

>>> dst_filename2 = "/tmp/xx_geotiff_copy2.tif"
>>> dst_ds = driver.CreateCopy( dst_filename2, src_ds, 0,
...    [ 'TILED=YES', 'COMPRESS=PACKBITS' ] )
>>> dst_filename3 = "/tmp/xx_geotiff_copy3.tif"
>>> dst_ds3 = driver.CreateCopy( dst_filename3, src_ds, 0,
...    [ 'TILED=YES', 'COMPRESS=PACKBITS' ] )

This example illustrates the use of tile storage instead of strip storage in the conversion of Tiff. Different software supports different storage methods. The compression method used is PACKBITS. The fourth parameter may have different options according to different formats, so it cannot be listed uniformly. Here you can see all the creation parameters of GTiff. The fifth and sixth parameters are to register a callback function that reflects the transformation process (such as drawing a progress bar).

Pixel storage order

TIFF format files are used more, and there is more to say about the pixel storage order. Pay more attention when using it. The CreateCopy()function called above has the keyword parameter TILED (here refers to the GeoTIFF file, some data formats cannot use this parameter.). Let’s take a look at how this parameter is not declared (running by default):

>>> driver.CreateCopy("/tmp/xx_geotiff_copy_a1.tif",src_ds,0)

The generated images are organized in bands, and GDAL sets the domain 284 of the TIFF file to 2 at the time of export.

If you add parameters, the following code:

>>> driver.CreateCopy("/tmp/xx_geotiff_copy_a2.tif",src_ds,0,["INTERLEAVE=PIXEL"])

The generated images are organized by pixels, and the specific storage method depends on the purpose of use. As for the way of organization, we can further understand the differences of BSQ format, BIL format and BIP format in remote sensing images.

Create Image Using Create Function

If you are not using an existing image file to create a new image, you need to use the Create method. In the data processing process, the Create method is the main method, which can output the virtual data set built in memory to the actual file. That is, the concept of raster data persistence, the in-memory data model (mainly two-dimensional array) is converted into a storage model. For GIS, in addition to the data itself, there are projections, metadata information, and so on.

The Createfunction is similar to the CreateCopy function, but it has several more parameters, including image size, band (channel) number, and pixel data type.

>>> import gdal
>>> driver = gdal.GetDriverByName( 'GTiff' )
>>> dst_filename = '/tmp/x_tmp.tif'
>>> dst_ds=driver.Create(dst_filename,512,512,1,gdal.GDT_Byte)

The above statement creates a single-band raster image in GeoTIFF format with a length and width of \(512\times 512\)pixels, and the data type is Byte. As an example, no other source data is used here, it creates an empty data set. Additional code is required if actual data is to be used.

>>> import numpy, osr
>>> dst_ds.SetGeoTransform([444720, 30, 0, 3751320, 0, -30 ])
>>> srs = osr.SpatialReference()
>>> srs.SetUTM( 11, 1 )
>>> srs.SetWellKnownGeogCS( 'NAD27' )
>>> dst_ds.SetProjection( srs.ExportToWkt() )
>>> raster = numpy.zeros( (512, 512) )
>>> dst_ds.GetRasterBand(1).WriteArray( raster )

The above example sets the spatial extent and coordinate system, and the size of \(512\times 512\) is all array data of 0.

Creat multiband images

Method of creating multi-band image

Most remote sensing images are multi-band. Each band records different spectral information. For image processing results, multi-band storage can also be used.

Here is a small example of building a 3-band GeoTiff. Multi-band images are created in a similar way.

>>> from osgeo import gdal
>>> import numpy
>>> dataset = gdal.Open("/gdata/geotiff_file.tif")
>>> width = dataset.RasterXSize
>>> height = dataset.RasterYSize
>>> datas = dataset.ReadAsArray(0,0,width,height)
>>> driver = gdal.GetDriverByName("GTiff")
>>> tods = driver.Create("/tmp/x_geotiff_file_3.tif",width,
...     height,3,options=["INTERLEAVE=PIXEL"])
>>> tods.WriteRaster(0,0,width,height,datas.tostring(),width,
...     height,band_list=[1,2,3])
0
>>> tods.FlushCache()

This is a very simple way to save remote sensing images (not including spatial information). In particular, pay attention to the options=["INTERLEAVE=PIXEL"] parameter in the Create() function. Note that this place, because the ``ReadAsArray()`` function is used when reading data, so when writing data to the tods variable, you need to use datas.tostring() to convert type of data. If you use the ReadRaster()function when reading data, you don’t need to convert it. It can be considered that ReadRaster() and WriteRaster() are corresponding/reciprocal; there is no function directly corresponding to ReadAsArray().

Also note the band_list parameter, which is determined by the number of bands. The parameter value needs to be judged based on the source data to prevent an abnormality.

Another problem is the order in which the data is stored, which is closely related to INTERLEAVE.

Bandwidth processing

It is also possible to read data from the band and then stitch it into a complete image. For splicing, you can use the numpy.concatenatefunction, which is a function provided by NumPy that can stitch multiple arrays at once.

>>> datas = []
>>> for i in range(3):
...     band = dataset.GetRasterBand(i+1)
...     data = band.ReadAsArray(0,0,width,height)
...     datas.append(numpy.reshape(data,(1,-1)))
>>> datas = numpy.concatenate(datas)

Then, output it as an image file:

>>> tods = driver.Create("/tmp/x_geotiff_file_4.tif",width,
...     height,3,options=["INTERLEAVE=PIXEL"])
>>> tods.WriteRaster(0,0,width,height,datas.tostring(),width,
...     height,band_list=[1,2,3])
0
>>> tods.FlushCache()

If you want each band to be input as a variable, you can loop to each band and then use the WriteRaster of the Band object instead of calling WriteRaster of Dataset.

Spatial Projection Processing in GDAL Writing Operation

Spatial reference and writing of projection information

When using GDAL to create data, it should be noted that the spatial parameter information of the image is processed separately. For example, importing a NAD27 spatial reference can be written as follows:

>>> import osr
>>> import gdal
>>> dataset = gdal.Open("/gdata/geotiff_file.tif")
>>> width = dataset.RasterXSize
>>> height = dataset.RasterYSize
>>> datas = dataset.ReadAsArray(0,0,width,height)
>>> driver = gdal.GetDriverByName("GTiff")
>>> tods = driver.Create("/tmp/x_geotiff_file_3.tif",width,
...     height,3,options=["INTERLEAVE=PIXEL"])
>>> tods.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )
0
>>> srs = osr.SpatialReference()
>>> srs.SetUTM( 11, 1 )
0
>>> srs.SetWellKnownGeogCS( 'NAD27' )
0
>>> tods.SetProjection( srs.ExportToWkt() )
0

So TIFF becomes a GeoTIFF file with a spatial reference system of NAD27, starting at (444720,3751320) and having a pixel size of 30. NAD27 is a common space reference code in the United States.

If the spatial reference used is troublesome, only six-parameter transformation can be defined.

Establishing the pyramids of images

The image pyramid is a set of images with different resolution from thin to thick, which is generated by the original image according to certain rules. At the bottom of the pyramid is the high resolution representation of the image, that is, the original image, while at the top is the approximation of the low resolution.

Generally, the amount of remote sensing image data is large. Pyramids are built for images, and these images can be displayed quickly. Avoiding the need to access the whole raster data set pyramid when displaying is a method of storing raster images in a copy mode that can reduce the resolution step by step.

In Python, image pyramids can be set up in different styles, among which Erdas Imagine is used more frequently.

>>> from osgeo import gdal
>>> gdal.SetConfigOption('HFA_USE_RRD', 'YES')

Establish image pyramids in the code:

>>> rds = gdal.Open("/gdata/lu75c.tif")
>>> rds.BuildOverviews(overviewlist=[2,4, 8,16,32,64,128])