Raster-Specific Operations

Raster-Specific Operations
Prev	Chapter 16. Temporal Tessellations	Next

The sections above state what a raster shares with the other tessellations. The operations below are the raster's own: a raster cell carries a value where an H3, QUADBIN or S2 cell carries only identity, so a trajectory read against a raster answers a tfloat rather than a temporal cell.

MobilityDB's raster sampling operator reads a raster band along a moving-point trajectory and returns a tfloat. A position outside the raster extent or on a nodata pixel carries no value. The sampling operators are double-valued whatever the pixel type of the band, so a pixel type belongs to the sampling surface when every value it can hold is exactly representable in a double precision number. That is what lets a band of any type be sampled into one tfloat, and a column holding tiles of several pixel types be read by a single query. A pixel type is also accepted under the name PostGIS raster gives it, so the band type an ST_BandPixelType call reports (8BUI, 16BSI, 32BF, …) passes into the tile constructors as it stands. PostGIS raster carries no 16-bit float and no 64-bit integer, so float16, int64 and uint64 take the spelling its grammar gives them, 16BF, 64BSI and 64BUI. Its 1BB, 2BUI and 4BUI types bound a pixel to less than a byte while PostGIS stores each of them a byte a pixel, so such a band is already the bytes of an unsigned 8-bit band and the constructors read it as uint8; the bound is the only thing the name adds, and a tile does not carry it. The name a tile reports is the one the RaQuet specification gives the type. The types uint64 and int64 hold values beyond that domain: a tile of either is stored and read back exactly, while sampling a pixel whose value no double precision number names raises an error rather than answering with a neighbouring value. This enables queries such as what elevation profile did each vehicle follow? or which trips entered a temperature threshold band?

A Raquet table keeps the QUADBIN cell of a tile in its block column as a 64-bit integer, beside the metadata row at block = 0, which holds no QUADBIN index. A query casts the cells of a trajectory to bigint to read that column, and a cell read from it to quadbin, the cast checking that it is a QUADBIN index. The typical query pattern joins a trajectory against a Raquet table in two steps:

-- Step 1: identify which tiles the trajectory overlaps
SELECT getValues(tquadbin(traj, 8));
-- Step 2: sample each matching tile, a 256 x 256 float32 band without nodata
SELECT rasterTileValueQuadbin(traj, band_1, 256, 256, block::quadbin, 'float32', 0,
  false) FROM elevation_raquet
WHERE block IN (SELECT unnest(getValues(tquadbin(traj, 8)))::bigint);

The restriction and predicate operators that close the section compose rasterValue with the standard MobilityDB temporal restriction and predicate functions. They accept an optional band argument (default 1) and call rasterValue exactly once per invocation.

A raster function takes the name of its PostGIS raster counterpart with the ST_ prefix removed, and the names and defaults of its arguments. One signature answers each operation, so a PostGIS form that only fixes an argument of another is written through that other form. A name the SQL standard reserves takes the type as a prefix, so ST_Value is rasterValue, and the Well-Known Binary functions take the names every MobilityDB type takes.

PostGIS raster	MobilityDB
`ST_Value`	`rasterValue`
`ST_Clip`	`clip`
`ST_Transform`	`transform`
`ST_Rescale`	`rescale`
`ST_Reclass`	`reclass`
`ST_SummaryStats`	`summaryStats`
`ST_DumpAsPolygons`	`dumpAsPolygons`
`ST_NumBands`, `ST_Width`, `ST_Height`, `ST_SRID`	`numBands`, `width`, `height`, `SRID`
`ST_UpperLeftX`, `ST_UpperLeftY`, `ST_ScaleX`, `ST_ScaleY`, `ST_SkewX`, `ST_SkewY`	`upperLeftX`, `upperLeftY`, `scaleX`, `scaleY`, `skewX`, `skewY`
`ST_BandPixelType`, `ST_BandNoDataValue`	`bandPixelType`, `bandNoDataValue`
`ST_AsBinary`, `ST_AsHexWKB`	`asBinary`, `asHexWKB`
`ST_RastFromWKB`, `ST_RastFromHexWKB`	`rasterFromBinary`, `rasterFromHexWKB`

asBinary and asHexWKB take a byte order where PostGIS takes outasin, and without one answer the bytes of ST_AsBinary and ST_AsHexWKB.

Read a raster band along a trajectory
```
rasterValue(tgeompoint,rast raster,band integer=1,exclude_nodata_value boolean=true,
  resample text='nearest') → tfloat
rasterValue(tgeompoint,path text,band integer=1) → tfloat
```
Reads band band along traj, taking the value of the pixel each position falls in. A trajectory that states nothing between its instants is read at its instants and answers an instant-set tfloat. A trajectory that moves between them passes over the pixels between them, so each segment is traversed pixel by pixel, and it answers a step tfloat holding each value from the instant the trip reaches the pixel holding it, a pixel whose corner the trip only clips included. Under an exclusive upper bound, a pixel the trajectory reaches only at its end, holding it for no time, carries no value into the result, so the reads of a trajectory cut at successive half-open periods merge into the read of the whole trajectory. A position outside the raster or over a nodata pixel carries no value and ends the run, so a trip that leaves and re-enters the raster answers one sequence per visit, and NULL when it never meets a pixel carrying data. A band whose pixels cannot be read, such as a band stored outside the database in a file that cannot be opened, raises an error rather than reading as a band without values. The form taking path reads a raster file on the server, in any format GDAL supports, through GDAL, where PostGIS allows a raster band stored outside the database to be read: postgis.enable_outdb_rasters must be on and postgis.gdal_enabled_drivers must enable the driver of the file. It answers what the form taking a raster answers over the same raster.
exclude_nodata_value and resample are those of the PostGIS ST_Value: with the first false a position over a nodata pixel keeps the nodata value where it is otherwise dropped, and the second reads the pixel under the position with nearest or bilinear, named without regard to case. A bilinear value varies within a pixel. Along a segment of a trajectory that moves between its instants it is, over each half-pixel square the segment crosses, a polynomial of degree two in time, so a trajectory read with bilinear answers a linear tfloat holding the value at the instants of the trajectory, where it crosses the pixel edges and the lines through the pixel centres, and at the extremum each of those pieces reaches, linear between them, as the product of two temporal floats holds its extremum. Its minimum and maximum are those of the value, and it is the value itself where the trajectory moves along a grid axis, over which the value is linear; elsewhere, between the instants it holds, it is the chord of the value. The value jumps across a pixel edge where a pixel read there holds nodata, which takes the value of the pixel the position falls in, and a new sequence then starts at the instant the trajectory crosses the edge.
```
-- Build a 3x3 synthetic elevation raster (SRID 4326, 1 degree pixels)
WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(3, 3, 0.0, 3.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, NULL::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4],
    [40.0::float4, 50.0::float4, 60.0::float4],
    [70.0::float4, 80.0::float4, 90.0::float4]]) AS r )
SELECT rasterValue(tgeompoint 'SRID=4326;[POINT(0.5 2.5)@2001-01-01,
  POINT(2.5 2.5)@2001-01-02, POINT(0.5 0.5)@2001-01-03]', r)::text FROM rast;
/* Interp=Step;[10@2001-01-01, 20@2001-01-01 06:00:00, 30@2001-01-01 18:00:00,
   50@2001-01-02 06:00:00, 70@2001-01-02 18:00:00, 70@2001-01-03] */
-- A bilinear read reaching its maximum between two pixel centres
WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(2, 2, 0.0, 2.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, NULL::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4], [20.0::float4, 10.0::float4]]) AS r )
SELECT round(rasterValue(tgeompoint 'SRID=4326;[POINT(0.5 1.5)@2001-01-01,
  POINT(1.5 1.1)@2001-01-03]', r, 1, true, 'bilinear'), 6)::text FROM rast;
/* [10@2001-01-01, 15@2001-01-02, 16.125@2001-01-02 18:00:00, 16@2001-01-03] */
```
Applied to a trajectory dataset with a terrain elevation raster:
```
-- Elevation profile per trip (requires a loaded elevation raster)
SELECT tripid, rasterValue(trip, elev) AS elevation FROM trips, elevation_raster;
-- Average elevation per trip
SELECT tripid, avgValue(rasterValue(trip, elev)) AS mean_elev FROM trips,
  elevation_raster;
-- Trips that crossed terrain above 200 m
SELECT tripid FROM trips, elevation_raster
WHERE atSpan(rasterValue(trip, elev), floatspan '[200, 9999]') IS NOT NULL;
```
A raster loaded as a table of tiles, as raster2pgsql -t loads it, is read along a trajectory by merging the reads of its tiles with mergeAgg. A position lies in a pixel of one tile, so the reads of the tiles hold disjoint times, and their merge answers the read of the raster the tiles cover, a trip crossing from one tile to the next included.
```
-- A 4x4 raster read along a trip whole, and as its four tiles of 2x2 pixels
WITH r AS (
  SELECT ST_SetValues(ST_AddBand(ST_MakeEmptyRaster(4, 4, 0.0, 4.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, NULL::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4, 40.0::float4],
    [50.0::float4, 60.0::float4, 70.0::float4, 80.0::float4],
    [90.0::float4, 100.0::float4, 110.0::float4, 120.0::float4],
    [130.0::float4, 140.0::float4, 150.0::float4, 160.0::float4]]) AS rast ),
t AS (
  SELECT tgeompoint 'SRID=4326;[POINT(0.5 3.5)@2001-01-01, POINT(3.5 0.5)@2001-01-02]'
    AS trip )
SELECT mergeAgg(rasterValue(trip, tile)) = (SELECT rasterValue(trip, rast) FROM t, r)
FROM t, r, ST_Tile(rast, 2, 2) AS tile;
-- true
```

Return a raster whose band holds the classes an expression maps its values to

reclass(raster,nband integer,reclassexpr text,pixeltype text,nodataval float=NULL)
  → raster

The expression follows the grammar of the PostGIS ST_Reclass: a comma-separated list of mappings, each a source range and a destination around a colon, a range being one value or two around a dash. A range written plainly includes its low end and excludes its high end, while ( excludes the low end and ] includes the high end, so 0-50 leaves 50 unmapped where [50-100] takes it. A pixel no range takes holds the nodata value of the resulting band, or zero when the band states none. That band has the pixel type pixeltype, named as PostGIS raster names it, such as 8BUI, and states no nodata value when nodataval is left out. The raster keeps its grid and its other bands, so the classes can be read along a trajectory with rasterValue.

-- Build a 3x3 synthetic elevation raster (SRID 4326, 1 degree pixels)
WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(3, 3, 0.0, 3.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, NULL::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4],
    [40.0::float4, 50.0::float4, 60.0::float4],
    [70.0::float4, 80.0::float4, 90.0::float4]]) AS r )
SELECT ST_DumpValues(reclass(r, 1, '0-50:1, 51-100:2', '8BUI', 0), 1) AS plain,
  ST_DumpValues(reclass(r, 1, '0-50:1, [50-100]:2', '8BUI', 0), 1) AS bracket
FROM rast;
-- {{1,1,1},{1,NULL,2},{2,2,2}} | {{1,1,1},{1,2,2},{2,2,2}}
-- Read the classes along a trajectory
WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(3, 3, 0.0, 3.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, NULL::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4],
    [40.0::float4, 50.0::float4, 60.0::float4],
    [70.0::float4, 80.0::float4, 90.0::float4]]) AS r )
SELECT rasterValue(tgeompoint 'SRID=4326;{POINT(0.5 2.5)@2001-01-01,
  POINT(1.5 1.5)@2001-01-02, POINT(2.5 0.5)@2001-01-03}',
  reclass(r, 1, '0-50:1, [50-100]:2', '8BUI', 0))::text FROM rast;
-- {1@2001-01-01, 2@2001-01-02, 2@2001-01-03}

Return a raster keeping the pixels of another that a geometry covers

clip(raster,geom geometry,crop boolean=true,touched boolean=false) → raster
clip(raster,nband integer[],geom geometry,nodataval float[]=NULL,crop boolean=true,
  touched boolean=false) → raster

Every band is read against a mask the geometry is burnt into, as the PostGIS ST_Clip reads it: a pixel the geometry does not cover answers the nodata value of its band, and the smallest value its pixel type holds when the band declares none. With crop the result carries the extent the raster and the geometry share, and otherwise the extent of the raster. The geometry is in the reference system of the raster; one in another system raises an error. With nband the result keeps the bands named, every band when it is NULL, nodataval states the nodata value of the pixels the geometry does not cover, one for every band or one per band, and with touched a pixel the geometry touches is kept where otherwise only one whose centre it covers is.

-- Build a 3x3 synthetic elevation raster (SRID 4326, 1 degree pixels)
WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(3, 3, 0.0, 3.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, -9999::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4],
    [40.0::float4, 50.0::float4, 60.0::float4],
    [70.0::float4, 80.0::float4, 90.0::float4]]) AS r )
SELECT ST_DumpValues(clip(r, geometry 'SRID=4326;POLYGON((0 1,2 1,2 3,0 3,0 1))'), 1)
  AS cropped,
  ST_DumpValues(clip(r, geometry 'SRID=4326;POLYGON((0 1,2 1,2 3,0 3,0 1))', false), 1)
  AS kept_extent
FROM rast;
-- {{10,20},{40,50}} | {{10,20,NULL},{40,50,NULL},{NULL,NULL,NULL}}

Return a raster stated in another spatial reference system

transform(raster,srid integer,algorithm text='NearestNeighbour',maxerr float=0.125,
  scalex float=0,scaley float=0) → raster
transform(raster,alignto raster,algorithm text='NearestNeighbour',
  maxerr float=0.125) → raster
transformPipeline(raster,pipeline text,srid integer=0,is_forward boolean=true,
  algorithm text='NearestNeighbour',maxerr float=0.125,scalex float=0,
  scaley float=0) → raster

Every band is carried into the target system and resampled onto the grid the reprojection implies, through the same GDAL warp as the PostGIS ST_Transform, so the result states the same coverage read through another system. The algorithm is one of NearestNeighbour, Bilinear, Cubic, CubicSpline, Lanczos, Max and Min, named without regard to case, and maxerr is the error in input pixels the warp may commit, 0 for an exact calculation. scalex and scaley state the pixel size of the result, 0 letting the warp derive it, and alignto takes the reference system, the pixel size, the grid origin and the skew of another raster, so the result lies on its grid.

The transformPipeline function carries the raster through a coordinate operation, stated as a PROJ string, a WKT or an authority code such as urn:ogc:def:coordinateOperation:EPSG::16031, instead of the one GDAL derives from the two systems, and reads it as the geometry transformPipeline does: longitude first and in degrees. The result states the reference system given by srid, 0 by default. By default the operation is applied forward; setting is_forward to false applies its inverse. The remaining arguments are those of transform, so an operation GDAL would derive answers as transform to the same system.

-- Build a 3x3 synthetic elevation raster (SRID 4326, 1 degree pixels)
WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(3, 3, 0.0, 3.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, -9999::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4],
    [40.0::float4, 50.0::float4, 60.0::float4],
    [70.0::float4, 80.0::float4, 90.0::float4]]) AS r )
SELECT SRID(transform(r, 3857)) AS srid, width(transform(r, 3857)) AS width,
  height(transform(r, 3857)) AS height
FROM rast;
-- 3857 | 3 | 3
-- Build a 3x3 synthetic elevation raster (SRID 4326, 1 degree pixels)
WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(3, 3, 0.0, 3.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, -9999::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4],
    [40.0::float4, 50.0::float4, 60.0::float4],
    [70.0::float4, 80.0::float4, 90.0::float4]]) AS r )
SELECT SRID(transformPipeline(r, '+proj=utm +zone=31 +ellps=WGS84', 32631)) AS srid,
  width(transformPipeline(r, '+proj=utm +zone=31 +ellps=WGS84', 32631)) AS width,
  height(transformPipeline(r, '+proj=utm +zone=31 +ellps=WGS84', 32631)) AS height
FROM rast;
-- 32631 | 3 | 3

Return a raster resampled to another pixel size

rescale(raster,scalex float,scaley float,algorithm text='NearestNeighbour',
  maxerr float=0.125) → raster

The result keeps the reference system and the upper left corner of the raster and states its coverage on a grid of the pixel size asked for, through the same GDAL warp as the PostGIS ST_Rescale, so a coarser scale reads fewer pixels over the same ground. The scales are in the units of the raster's reference system, and the algorithm and the error are those of transform.

-- Build a 3x3 synthetic elevation raster (SRID 4326, 1 degree pixels)
WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(3, 3, 0.0, 3.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, -9999::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4],
    [40.0::float4, 50.0::float4, 60.0::float4],
    [70.0::float4, 80.0::float4, 90.0::float4]]) AS r )
SELECT width(rescale(r, 0.5, -0.5)) AS width,
  height(rescale(r, 0.5, -0.5)) AS height,
  scaleX(rescale(r, 0.5, -0.5)) AS scalex
FROM rast;
-- 6 | 6 | 0.5

Return the statistics of the pixels of a raster band

summaryStats(raster,nband integer=1,exclude_nodata_value boolean=true) → summarystats

The band is read once and answers the PostGIS summarystats record that ST_SummaryStats answers: the number of pixels counted, and the sum, mean, standard deviation, minimum and maximum of their values. Every pixel is read, so the answer is exact rather than sampled; the pixels the band states as nodata are left out unless exclude_nodata_value is false, and a band with no pixel to count answers a count of zero and no statistic.

-- Build a 3x3 synthetic elevation raster (SRID 4326, 1 degree pixels)
WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(3, 3, 0.0, 3.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, -9999::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4],
    [40.0::float4, 50.0::float4, 60.0::float4],
    [70.0::float4, 80.0::float4, 90.0::float4]]) AS r )
SELECT (summaryStats(r)).* FROM rast;
-- 9 | 450 | 50 | 25.81988897471611 | 10 | 90

Return the polygons of a raster band, one for each group of pixels sharing a value

dumpAsPolygons(raster,band integer=1,exclude_nodata_value boolean=true) → {geomval}

A band states a value per pixel, and this states the same band as the regions those values cover, answering the PostGIS geomval records that ST_DumpAsPolygons answers. Each polygon carries the reference system of the raster, where ST_DumpAsPolygons states none, so it can be compared with the trajectories the coverage is read along. The pixels the band states as nodata are left out unless exclude_nodata_value is false, a band whose every pixel is nodata covers nothing and answers no row, and a band the raster does not have raises an error.

-- Classify a 3x3 synthetic elevation raster, then read each class as a region
WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(3, 3, 0.0, 3.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, -9999::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4],
    [40.0::float4, 50.0::float4, 60.0::float4],
    [70.0::float4, 80.0::float4, 90.0::float4]]) AS r )
SELECT val, ST_AsText(geom)
FROM rast, dumpAsPolygons(reclass(r, 1, '0-50:1, [50-100]:2', '8BUI', 0))
ORDER BY val;
-- 1 | POLYGON((0 3,0 1,1 1,1 2,3 2,3 3,0 3))
-- 2 | POLYGON((1 2,1 1,0 1,0 0,3 0,3 2,1 2))

Sample a raquet raster tile along a trajectory
```
rasterTileValue(tgeompoint,rast raquet) → tfloat
```
Reads the tile along traj (SRID 4326) as rasterValue reads a raster, returning the pixel values as a tfloat. This is the typed equivalent of rasterTileValueQuadbin: the dimensions, QUADBIN cell, pixel type and nodata handling are taken from the raquet value instead of from separate arguments. Returns NULL when the trip never meets a pixel of the tile carrying data.
```
-- Sample pre-packaged raquet tiles along ship trajectories
SELECT t.tripid, rasterTileValue(t.trip, r.tile)
FROM trips t JOIN elevation_tiles r
  ON r.quadbin IN (SELECT unnest(getValues(tquadbin(t.trip, 8))));
```
Sample an array of raquet raster tiles along a trajectory
```
rasterTileValue(tgeompoint,rast raquet[]) → tfloat
```
Reads every tile of rast along traj (SRID 4326) and returns the pixel values as a single tfloat. A trajectory that leaves one tile is covered by several of them, each answering the part of the trip over it. Tiles of one zoom level partition the plane, so they answer disjoint times. Tiles of different zoom levels overlap, and where two of them answer the same time the value of the tile of higher zoom is kept, that being the one carrying the finer resolution. Returns NULL when the trip never meets a pixel of a tile carrying data.
```
-- Sample a trajectory across every tile it crosses, in one value
SELECT t.tripid, rasterTileValue(t.trip, array_agg(r.tile)) FROM trips t
JOIN elevation_tiles r ON r.quadbin IN (SELECT unnest(getValues(tquadbin(t.trip, 8))))
GROUP BY t.tripid, t.trip;
```

Sample a Raquet raster chip along a trajectory using a QUADBIN cell for georeferencing

rasterTileValueQuadbin(tgeompoint,pixels bytea,width integer,height integer,
  quadbin,pixtype text,nodata float,has_nodata boolean) → tfloat

Reads a row-major single-band pixel array along traj (SRID 4326), as rasterValue reads a raster. The tile georeferencing (bounding box, pixel-to-coordinate mapping) is derived from quadbin alone via the Web-Mercator slippy-tile transform. The pixtype argument must be one of uint8, int8, uint16, int16, uint32, int32, uint64, int64, float16, float32, or float64. When has_nodata is true, a pixel equal to nodata carries no value. Returns NULL when the trip never meets a pixel carrying data.

-- Sample the 256 x 256 float32 elevation tiles of a Raquet table, stating no
-- nodata value, along every ship trajectory
SELECT t.tripid, rasterTileValueQuadbin(t.trip, r.band_1, 256, 256,
  r.block::quadbin, 'float32', 0, false) FROM trips t
JOIN elevation_raquet r
  ON r.block IN (SELECT unnest(getValues(tquadbin(t.trip, 8)))::bigint);
-- Average sea-surface temperature per trajectory
SELECT t.tripid, avgValue(rasterTileValueQuadbin(t.trip, 
  r.band_1, 256, 256, r.block::quadbin, 'int16', -9999, true)) AS mean_sst FROM trips t
JOIN sst_raquet r ON r.block IN (SELECT unnest(getValues(tquadbin(t.trip, 6)))::bigint);

Return the part of a trajectory where the raster value it reads falls inside a float range

atRasterValue(tgeompoint,rast raster,vspan floatspan,band integer=1) → tgeompoint
atRasterValue(tgeompoint,path text,vspan floatspan,band integer=1) → tgeompoint

Reads rasterValue along traj and returns the sub-trajectory restricted to the times when the pixel value lies within vspan. Returns NULL when no time satisfies the condition. The form taking path reads a raster file on the server, as the corresponding form of rasterValue does.

WITH rast AS (
  SELECT ST_SetValues(ST_AddBand( ST_MakeEmptyRaster(3, 3, 0.0, 3.0, 1.0, -1.0, 0.0, 0.0,
    4326), '32BF'::text, 0.0::float8, NULL::float8), 1, 1, 1,
    ARRAY[[10.0::float4, 20.0::float4, 30.0::float4],
    [40.0::float4, 50.0::float4, 60.0::float4],
    [70.0::float4, 80.0::float4, 90.0::float4]]) AS r )
SELECT asText(atRasterValue(tgeompoint 'SRID=4326;
  [POINT(0.5 2.5)@2001-01-01, POINT(1.5 1.5)@2001-01-02, POINT(0.5 0.5)@2001-01-03]', r,
  floatspan '[40, 90]')) FROM rast;
/* {[POINT(1 2)@2001-01-01 12:00:00, POINT(1.5 1.5)@2001-01-02,
   POINT(0.5 0.5)@2001-01-03]} */

Return the part of a trajectory where the raster value it reads falls outside a float range
```
minusRasterValue(tgeompoint,rast raster,vspan floatspan,band integer=1) → tgeompoint
minusRasterValue(tgeompoint,path text,vspan floatspan,band integer=1) → tgeompoint
```
The complement of atRasterValue: returns the sub-trajectory restricted to the times when the pixel value lies outside vspan. Returns NULL when the pixel value lies within vspan at every time. The form taking path reads a raster file on the server, as the corresponding form of rasterValue does.
```
-- Keep only the instant whose raster value (10) is below 40
SELECT asText(minusRasterValue(traj, elev_raster, floatspan '[40, 9999]')) FROM trips,
  elevation_raster;
```

Return true if the trajectory ever reads a raster pixel value inside a float range

eRasterValue(tgeompoint,rast raster,vspan floatspan,band integer=1) → boolean
eRasterValue(tgeompoint,path text,vspan floatspan,band integer=1) → boolean

Returns true when traj ever reads a pixel value within vspan, false otherwise. Returns NULL only when the input is NULL. The form taking path reads a raster file on the server, as the corresponding form of rasterValue does.

-- Trips that ever crossed terrain above 200 m
SELECT tripid FROM trips,
  elevation_raster WHERE eRasterValue(trip, elev_raster, floatspan '[200, 9999]');

Return true if every pixel value the trajectory reads over the raster lies inside a float range
```
aRasterValue(tgeompoint,rast raster,vspan floatspan,band integer=1) → boolean
aRasterValue(tgeompoint,path text,vspan floatspan,band integer=1) → boolean
```
Returns true when every pixel value traj reads over the raster lies within vspan, false when at least one does not. Returns NULL only when the input is NULL. The form taking path reads a raster file on the server, as the corresponding form of rasterValue does.
```
-- Trips that stayed entirely below 100 m elevation
SELECT tripid FROM trips,
  elevation_raster WHERE aRasterValue(trip, elev_raster, floatspan '[0, 100]');
```

To use the operators described here, MobilityDB must be built with -DQUADBIN=ON and -DRASTER=ON, since a Raquet tile is identified by a quadbin cell. On such a build the generated mobilitydb.control declares requires = 'postgis, postgis_raster', so a single CASCADE creates the full stack:

CREATE EXTENSION mobilitydb CASCADE;
-- NOTICE:  installing required extension "postgis"
-- NOTICE:  installing required extension "postgis_raster"

PostGIS raster is part of the standard postgresql-N-postgis-3 package on Debian/Ubuntu; no extra package is required. Pass -DQUADBIN=ON and -DRASTER=ON to CMake:

cmake -DQUADBIN=ON -DRASTER=ON ..
make -j$(nproc)
sudo make install

The CI builds, .github/workflows/pgversion.yml among them, pass -DALL=ON, which includes raster alongside the other optional families.

Prev	Up	Next
S2-Specific Operations	Home	Chapter 17. Temporal Point Clouds