Raster Tutorial

ISCE delegates raster data handling to the Geospatial Data Abstraction Library (GDAL) at the C++ level via isce::io::Raster. Even the Python interface to isce::io::Raster is designed to work in combination with GDAL’s Python bindings. We will walk through some of the examples here.

A concise introduction to GDAL’s Python bindings can be found here. Familiarizing yourself with this will help navigate the following tutorial more easily. Note that pyRaster is a thin wrapper over GDAL Dataset and is only needed when you want to pass on the created object to processing modules in ISCE. For simple I/O at the python level, one should directly use the GDAL interface.

  1. Open an existing file in Read-Only Mode

  2. Open an existing file in Update Mode

  3. Creating a raster

  4. Creating a raster using numpy array

  5. Creating a raster using h5py

Open an existing file in Read-Only Mode

from isce3.io import raster

#Create object
image = raster(filename="inputfilename")
print('Dims: {0}P x {1}L'.format(raster.width, raster.length))

#image is ready to be passed on to ISCE processing modules

Open an existing file in Update Mode

from isce3.io import raster
from osgeo import gdal

#Create object
ds = gdal.Open(inputfilename, gdal.GA_Update)
image = raster(dataset=ds)
for ii in range(image.numBands):
   print('Band {0} is of type {1}'.format(ii+1, gdal.GetDataTypeName(image.getDatatype(ii+1))))

 #image is ready to be passed on to ISCE processing modules

Creating a raster

It is really easy to create raster with GDAL and then pass it to ISCE to work with.

from isce3.io import raster
from osgeo import gdal

#Create GDAL raster
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create("antarctica.tif", 1000, 1500, 2, gdal.GDT_Float32);
ds.SetGeoTransform([-1.0e6, 1.0e3, 0., 1.5e6, 0., -1.0e3])

#Wrap it with pyRaster
image = raster(dataset=ds)

#Set projection code. Can do this with GDAL+osr as well before creating raster.
image.EPSG = 3031

#image is ready to be passed on to ISCE processing modules

Creating a raster using numpy array

You can also create GDAL datasets out of numpy arrays and pass it to ISCE to work with.

from isce3.io import raster
from osgeo import gdal_array
import numpy as np

#Create numpy array
arr = np.ones((1500,1000), dtype=np.complex64)

#Dress numpy array with gdalarray
ds = gdal_array.OpenArray(arr)

#Pass gdal dataset to pyRaster
image = raster(dataset=ds)

#image is ready to be passed on to ISCE processing modules

Creating a raster using h5py

You can also create ISCE Rasters out of h5py datasets. Note that for read only operations, you can use GDAL’s <a href=”https://gdal.org/drivers/raster/hdf5.html”>HDF5 driver</a> as well and set up rasters as shown above.

from isce3.io import raster
import h5py

#Create HDF5 file
fid = h5py.File('example.h5', 'w')

#Create group
grp = fid.create_group('level1/level2')

#Create dataset
#All complicated creation options can be controlled via h5py if needed
dset = grp.create_dataset("data", shape=(100,150), dtype='f4')

#Wrap it with ISCE raster
image = raster(h5=dset)

#image is ready to be passed on to ISCE processing modules