Introduction

This section outlines the general usage of the Python bindings for a Raster.

Reading and Writing of Rasters uses Numpy arrays so numpy must be installed for the version of Python being used.

All the examples assume the following imports has been done.

from caris.coverage import *
import numpy

Checking the type of a dataset

You can use caris.coverage.identify() to determine the type of dataset a file contains.

from caris.coverage import *

file_path = r'C:\bathymetry.csar'
type = identify(file_path)

if type == DatasetType.RASTER:
    raster = Raster(file_path)

elif type == DatasetType.CLOUD:
    cloud = Cloud(file_path)

elif type == DatasetType.VRS:
    vrs = VRS(file_path)

Opening a Raster

The following will open a raster for reading.

raster = Raster('source_raster.csar')

One can specify the named parameter.

raster = Raster(filename='source_raster.csar')

One can also supply a URI

raster = Raster(uri='file:///source_raster.csar')

Alternatively, the caris.open() function can be used:

from caris.coverage import *
import caris

# open file read only
raster = caris.open(file_name='raster.csar')

# open file read write
raster = caris.open(file_name='raster.csar', open_mode=caris.OpenMode.READ_WRITE)

Opening from BDB Server

To open a database surface, a URI must be given formatted as:

bdb://username:password@hostname/database/boid

Note

The database name is case sensitive.

Boids (object IDs for database objects) can be found by via caris.bathy.db.DatabaseFeature.id.

In BASE Editor the boid can be found by selecting a surfac, the boid will be shown in the Selection window. The URI for opened rasters or clouds is shown in the Properties window’s “Surface Name” field.

raster = caris.coverage.Raster(uri='bdb://dba:sql@example.com/MyDB/9245b052-634c-11e7-8dcf-1866da47ee87')

Listing Bands

The band information is found in the Raster.band_info property. It is simply a dictionary where the key is the band name and the value is an instance of BandInfo. The band information contains properties such as its type, category, minimum, maximum, no-data-value, etc. See BandInfo for more info.

raster = Raster('source_raster.csar')

for band_name in raster.band_info:
    print(str(raster.band_info[band_name]))

Reading from a Raster

Reading is done via the read() method.

The following will read the entire band named Depth. Note that this is not a good idea for large rasters. Reading in smaller blocks would be advisable.

raster = Raster('source_raster.csar')

data = raster.read(band_name = 'Depth', area = ((0,0),raster.dims), level = raster.highest_level)

The level is optional. It defaults to the highest level. This is equivalent code.

raster = Raster('source_raster.csar')

data = raster.read(band_name = 'Depth', area = ((0,0),raster.dims))

The levels will be numerically lower for the coarser levels down to raster.lowest_level inclusive.

Note that in most cases, using named parameters is an optional feature. One may specify the arguments in order without the names.

raster = Raster('source_raster.csar')

data = raster.read('Depth', ((0,0),raster.dims))

Reading from a Raster into a Numpy array

Rasters can also be read into a Numpy array via the read_np_array.

raster = Raster('source_raster.csar')

data = raster.read_np_array(band_name = 'Depth', area = ((0,0),raster.dims), level = raster.highest_level)

for x in range(0, data.shape[0] - 1):
    for y in range(0, data.shape[1] - 1):
        depth = data[x][y]

Example counting no-data-values

raster = Raster('source_raster.csar')

area = ((0,0), raster.dims)
ndvCount = sum(1 for x in raster.read('Depth', area) if x == raster.band_info['Depth'].ndv)
print('no-data-value count: ' + str(ndvCount))

# Or using numpy array
data = raster.read_np_array('Depth', area)
ndvCount = sum(1 for x in data.reshape(data.shape[0] * data.shape[1]) if x == raster.band_info['Depth'].ndv)

Examples converting between grid and geographic coordinates

Grid coordinates go from (0,0) to Raster.dims. They are in pixels. Suppose you wanted to know what the coordinate (100,100) is in geographic coordinates. Use the convertGridToGeo() method on the Raster’s transform.

raster = Raster('source_raster.csar')

geoPoint = raster.transform.convertGridToGeo((100,100), raster.highest_level)

Converting a grid area is similar.

raster = Raster('source_raster.csar')

gridArea = ((100,100),(200,200))
geoArea = raster.transform.convertGridToGeo(gridArea, raster.highest_level)

To convert in the other direction, use the convertGeoToGrid() method.

raster = Raster('source_raster.csar')

gridPoint = raster.transform.convertGeoToGrid((417760.1, 5579348.0), raster.highest_level)

raster = Raster('source_raster.csar')

geoArea = ((417760.1, 5579348.0),(417770.1, 5579358.0))
gridArea = raster.transform.convertGeoToGrid(geoArea, raster.highest_level)

Examples converting areas to a different level

If you have an area defined at the highest resolution, but you now want to read from a coarser level that represents the same geographic area, you can simply convert your coordinates before reading.

raster = Raster('source_raster.csar')

gridArea = ((100,100),(200,200))
gridArea2 = raster.transform.convertGridBox(raster.highest_level, raster.highest_level - 1, gridArea)
data = raster.read('Depth', gridArea2, raster.highest_level - 1)

Example expanding a geographic area to pixel boundary

Geographic coordinates don’t necessarilly match up to pixel coordinates. It is often desirable to expand the geographic extents to fully encompass all of the underlying pixels. The following example demonstrates how to do such a read.

raster = Raster('source_raster.csar')

geoArea = ((417760.1, 5579348.0),(417770.1, 5579358.0))

# Expand extents before converting to grid coordinates.
geoArea = raster.transform.expandGeo(geoArea, raster.highest_level)

# Now convert to grid space
gridArea = raster.transform.convertGeoToGrid(geoArea, raster.highest_level)

# Finally, read our data
data = raster.read('Depth', gridArea)

Copying a Raster

There are many ways to copy a Raster. The simplest way is to call the create_copy() member function.

raster = Raster('source_raster.csar')
outputFilename = 'dest_raster.csar'
raster.create_copy(outputFileName)

Another way of doing it is opening a Raster and copying all the required info into an Options instance.

raster = Raster('source_raster.csar')
outputFilename = 'dest_raster.csar'

# Setup options
opts = Options()
opts.open_type = OpenType.WRITE
opts.band_info = raster.band_info
opts.dims = raster.dims
opts.wkt_cosys = raster.wkt_cosys
opts.iso19139_xml = raster.iso19139_xml
opts.raster2geo_matrix = raster.transform.matrix
opts.pixel_type = raster.transform.pixel_type
opts.raster_read_callback_func = raster.read

# Create Raster
raster2 = Raster(outputFilename, options = opts)

Converting between BandInfo and a dictionary

A BandInfo class is not a dictionary and is not modifiable. However, one can convert to a dictionary, apply the changes you want and then convert it back. Here is how to convert to a dictionary.

raster = Raster('source_raster.csar')
bi = raster.band_info['Depth']

# Do the conversion
dbi = dict((x, getattr(bi, x)) for x in bi)

Converting back is much easier. Simply pass the dictionary to the BandInfo‘s constructor.

raster = Raster('source_raster.csar')
bi = raster.band_info['Depth']

# Do the conversion
dbi = dict((x, getattr(bi, x)) for x in bi)

# Convert back
bandInfo = BandInfo(dbi)

Example modifying Level Policies

A Level Policy is how to build the other pyramid levels of the raster. Each band can have its own level policy. Here is an example where we copy a Raster with known bands and we modify each band’s level policy.

raster = Raster('source_raster.csar')
outputFilename = 'dest_raster.csar'

# Setup options
opts = Options()
opts.open_type = OpenType.WRITE
opts.band_info = raster.band_info

# Change level policies
bi = opts.band_info['Deep']
dbi = dict((x, getattr(bi, x)) for x in bi)
dbi['level_policy'] = LevelPolicy.MIN
opts.band_info['Deep'] = BandInfo(dbi)

bi = opts.band_info['Density']
dbi = dict((x, getattr(bi, x)) for x in bi)
dbi['level_policy'] = LevelPolicy.FOLLOW
dbi['follow_band_name'] = 'Deep'
opts.band_info['Density'] = BandInfo(dbi)

bi = opts.band_info['Mean']
dbi = dict((x, getattr(bi, x)) for x in bi)
dbi['level_policy'] = LevelPolicy.FOLLOW
dbi['follow_band_name'] = 'Deep'
opts.band_info['Mean'] = BandInfo(dbi)

bi = opts.band_info['Shoal']
dbi = dict((x, getattr(bi, x)) for x in bi)
dbi['level_policy'] = LevelPolicy.MAX
opts.band_info['Shoal'] = BandInfo(dbi)

bi = opts.band_info['Weight']
dbi = dict((x, getattr(bi, x)) for x in bi)
dbi['level_policy'] = LevelPolicy.FOLLOW
dbi['follow_band_name'] = 'Shoal'
opts.band_info['Weight'] = BandInfo(dbi)

# Setup remaining options for copying
opts.dims = raster.dims
opts.wkt_cosys = raster.wkt_cosys
opts.iso19139_xml = raster.iso19139_xml
opts.raster2geo_matrix = raster.transform.matrix
opts.pixel_type = raster.transform.pixel_type
opts.raster_read_callback_func = raster.read

# Create raster
raster2 = Raster(outputFilename, options = opts)

Copying a Raster with Lambda

When creating a Raster, we must always copy from another source. This need not be from a Raster. In this case, we will use a lambda function. It does read from a source raster, but can be redirected to whatever source the user wishes.

def proxy_read_function(raster, offset, *args, **kwargs):
    in_box = args[1]
    level = raster.highest_level
    if (len(args) >= 3):
        level = args[2]
    if ('level' in kwargs):
        level = kwargs['level']

    box = ((in_box[0][0] + offset[0], in_box[0][1] + offset[1]), (in_box[1][0] + offset[0], in_box[1][1] + offset[1]))
    return raster.read(args[0], box, level)

def create_read_function(raster, offset):
    func = lambda args, kwargs, raster=raster, offset=offset: proxy_read_function(raster, offset, args, kwargs)
    return func

def create_copy_with_lambda():
    raster = Raster('source_raster.csar')
    outputFilename = 'dest_raster.csar'

    # Setup options
    opts = Options()
    opts.open_type = OpenType.WRITE
    opts.band_info = raster.band_info

    opts.dims = raster.dims
    opts.wkt_cosys = raster.wkt_cosys
    opts.iso19139_xml = raster.iso19139_xml
    opts.raster2geo_matrix = raster.transform.matrix
    opts.pixel_type = raster.transform.pixel_type

    # Use lambda function for reading from source.
    opts.raster_read_callback_func = create_read_function(raster, (0,0))

    # Create raster
    raster2 = Raster(outputFilename, options = opts)

Writing to a Raster

A Raster can be updated via Raster.write

import numpy

opts = Options()
opts.open_type = OpenType.WRITE

raster = Raster('source_raster.csar', options=opts)
area_to_write = ((0, 0), raster.dims)
band_name = 'Depth'
band_dtype = raster.band_info[band_name].numpy_dtype

# create some random data
num_values = raster.dims[0] * raster.dims[1]
array = numpy.arange(num_values, dtype=band_dtype).reshape(raster.dims[1], raster.dims[0])
# note raster.dims is columns x rows while numpy shapes are rows x columns

# write the data
raster.write(band_name, area_to_write, array)

Creating a Raster

A raster can be created from scratch via caris.coverage.create_raster().

uri = 'raster.csar'
origin = [0.0, 0.0]
resolution = [5.0, 5.0]
dimensions = [4, 4]

crs = """
PROJCS["WGS 84 / World Mercator",
GEOGCS["WGS 84",
        DATUM["World Geodetic System 1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
AUTHORITY["EPSG","7030"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
UNIT["degree (supplier to define representation)",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"],
EXTENSION["tx_authority","WG84"]],
PROJECTION["Mercator_1SP",
AUTHORITY["EPSG","19883"]],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","3395"]]
"""

band_info = BandInfo(name="Depth",
                     type = DataType.FLOAT32,
                     tuple_length = 1,
                     direction = Direction.DEPTH,
                     units = 'm',
                     category = Category.SCALAR,
                     level_policy = LevelPolicy.MAX)

bands = [band_info]
raster = create_raster(uri, crs, origin,
                       resolution, dimensions,
                       bands)

# write data
band_dtype = raster.band_info['Depth'].numpy_dtype
area = ((0, 0),(4, 4))

array = numpy.arange(16, dtype=band_dtype).reshape(4, 4)
raster.write("Depth", area, array)

Adding and removing bands

Bands can be added and removed with Raster.add_band and Raster.remove_band

import numpy
from caris.coverage import *

opts = Options()
opts.open_type = OpenType.WRITE

raster = Raster('raster.csar', options=opts)

band_name = "Test"

band_info = BandInfo(name=band_name,
type = DataType.FLOAT32,
tuple_length = 1,
direction = Direction.DEPTH,
units = 'm',
category = Category.SCALAR,
level_policy = LevelPolicy.MAX)

# add the band
raster.add_band(band_info)

# prepare some data to write
area_to_write = ((0, 0), raster.dims)
band_dtype = raster.band_info[band_name].numpy_dtype

# create some random data
num_values = raster.dims[0] * raster.dims[1]
array = numpy.arange(num_values, dtype=band_dtype).reshape(raster.dims[1], raster.dims[0])
# note raster.dims is columns x rows while numpy shapes are rows x columns

# write the data
raster.write(band_name, area_to_write, array)

# remove the band
raster.remove_band(band_name)

Updating bounding polygon

A bounding polygon can be automatically generated for a Raster. Any valid polygon can be set to be used as the bounding polygon. Used together, a bounding polygon can be replaced with the automatically generated bounding polygon.

opts = Options()
opts.open_type = OpenType.WRITE

raster = Raster(file_path, options=opts)
new_bounding_polygon = generate_polygon(raster)
raster.bounding_polygon = new_bounding_polygon