Module 5, Basic analyses with grid data

learning objective: learning how to work with grid domains, convert between vector and grid representations, perform neighbourhood operations, and navigate grid coordinates

introduction

Grid data is fundamentally different from vector data. In vector data, each feature has an explicit geometry attribute that records its location. In grid data, there is no geometry attribute at all: the spatial position of each cell is implicit in the domain unit itself. The grid domain defines a two-dimensional matrix of cells, and a projection links those cell positions to real-world coordinates.

This has two important consequences. First, calculations on grids are typically very fast because the regular structure allows specialised algorithms (such as FFT-based convolutions). Second, linking grid data to vector data always involves an explicit conversion step: you need to either rasterise your vector data onto the grid, or assign grid cells to vector features.

This module covers the most important operations for grid analysis in GeoDMS:

  • understanding and defining a grid domain
  • converting vector data to a grid (poly2grid)
  • grid coordinate functions (pointcol, pointrow)
  • neighbourhood operations (potential, proximity, diversity)
  • connected region analysis (district)
  • grid-to-vector aggregation

All examples use a 100-metre land-use grid (bbg2015_100m) together with the COROP region polygons from Module 2A.

the grid domain

A grid domain is a two-dimensional domain unit. Its value type is a point type such as spoint (signed 16-bit integer coordinates), ipoint (signed 32-bit), wpoint (unsigned 16-bit) or upoint (unsigned 32-bit). The most common choice for Dutch national grids is spoint for grids up to roughly 32,000 x 32,000 cells, or ipoint for larger ones.

The spatial extent of the grid — which real-world area it covers and at what resolution — is defined by the range property and the projection stored in the source file (for GeoTiff) or configured explicitly.

When you read a GeoTiff in Module 2A, the grid domain was defined implicitly by the file:

unit<spoint> land_use
:   StorageName     = "%ProjDir%/Data/bbg2015_100m_10k.tif"
,   StorageType     = "gdal.grid"
,   StorageReadOnly = "True"
,   DialogData      = "rdc_base"
{
   attribute<uint8> GridData;
}

You can also define a grid domain explicitly, which is useful when you want to create a new grid at a specific resolution and extent, or combine multiple rasters on a shared domain:

// Define a 100m grid covering the Netherlands in RD New
unit<spoint> grid_100m := range(spoint, point(-50s, 0s), point(625s, 300s))
{
   // The projection: origin in RD New, cell size 100m
   // (configured via the coordinate system unit, see the Grid-Domain page on the GeoDMS wiki)
}

Reading tip: for a full explanation of how to configure grid domains and projections explicitly, see grid domain

poly2grid: rasterising vector data

poly2grid(polygon_attribute, grid_domain) converts polygon vector data to a grid. For each cell in the grid domain, it returns the index of the polygon that covers that cell. The result is a relation from the grid domain to the polygon domain, and it works exactly like any other relation in GeoDMS: you can use it for lookups and aggregations.

// Rasterise COROP regions onto the 100m land-use grid
attribute<CBS_COROP> corop_grid (land_use) :=
    poly2grid(CBS_COROP/geometry, land_use);

The result corop_grid has the land_use grid as its domain and CBS_COROP as its values unit. Each cell gets the index of the COROP region it falls in, or null if no polygon covers it.

Once you have this relation, you can use it just like any other relation:

// Bring the COROP name to each grid cell
attribute<string> corop_name_grid (land_use) := CBS_COROP/statnaam[corop_grid];

// Aggregate grid cell counts per COROP region
attribute<uint32> n_cells_per_corop (CBS_COROP) := sum(const(1[uint32], land_use), corop_grid);

// Aggregate a grid value: the mean land-use code per COROP
// (illustrative only; land-use codes are categorical)
attribute<float32> mean_landuse_per_corop (CBS_COROP) :=
    mean(float32(land_use/GridData), corop_grid);

Performance note: poly2grid speed is primarily determined by the number of vertices in the polygon, not the number of cells. Simplifying complex polygons before rasterising can significantly improve performance. See geos_simplify_multi_polygon in Module 4.

poly2grid assigns at most one polygon per cell. If polygons overlap, the first one found is used. For a weighted cross-table that records the fractional overlap of each polygon with each cell, use poly2allgrids instead.

grid coordinate functions

Grid cells are addressed by their position in the two-dimensional domain. There are two ways to query this position, and they give fundamentally different results.

pointcol and pointrow: local grid position

pointcol(grid_attribute) and pointrow(grid_attribute) return the column and row number of each cell within the grid’s own local coordinate system. These are integer values starting from zero at the top-left corner of the grid.

attribute<int16> col (land_use) := pointcol(id(land_use));
attribute<int16> row (land_use) := pointrow(id(land_use));

The function id(domain) generates an attribute of the domain over itself, i.e. the identity: each element maps to its own index position. For grid domains, this is the cell’s 2D index.

Local coordinates are useful for position-relative calculations: finding neighbours, defining kernels, or encoding a cell’s position as a row/column pair for export.

geographic coordinates

get_x and get_y are the functions for retrieving the actual geographic (RD New) coordinates of each cell’s centre point.

neighbourhood operations

Neighbourhood operations compute, for each cell, a summary value based on the cells within a defined neighbourhood around it. In other GIS software these are sometimes called focal statistics or moving window operations.

the kernel

The neighbourhood is defined by a kernel: a small grid domain (not geographically projected) where each cell holds a weight. The kernel is centred on each input cell in turn, the input cell values are multiplied by the kernel weights, and the results are summed.

A simple example: a 3x3 uniform kernel (all weights 1) computes the sum of all 8 neighbours plus the cell itself.

The kernel is always defined as a separate, non-geographic grid domain with an spoint value type. The range defines the kernel’s extent, symmetrically around the origin (0, 0):

// A circular kernel with radius 5 cells (in 100m grid units: 500m radius)
unit<spoint> kernel_500m := range(spoint, point(-5s, -5s), point(6s, 6s))
{
   // dist2 gives the squared distance from each cell to the kernel centre
   attribute<uint32> dist2_val := dist2(point(0s, 0s, kernel_500m), uint32);

   // Distance-decay weight: 1.0 within 500m, 0 outside
   attribute<float32> weight := dist2_val <= 25 ? 1.0f : 0.0f;
}

The dist2 function (squared distance) is used here rather than sqrt(dist2) because it avoids the square root computation, which is expensive for large kernels. For a circular neighbourhood, the condition dist2 <= r^2 is equivalent to dist <= r.

potential: weighted sum over neighbourhood

potential(grid_attribute, kernel_weights) computes a neighbourhood sum: for each cell, the values of all cells within the kernel’s reach are multiplied by the corresponding kernel weight and summed. The result is a float32 grid.

// Count the number of inhabited cells (land-use class 1) within 500m of each cell
attribute<bool>    is_residential (land_use) := land_use/GridData == 1b;

attribute<float32> residential_potential (land_use) :=
    potential(float32(is_residential), kernel_500m/weight);

The result gives, for each 100m cell, the sum of weights of all residential cells within 500m. With a uniform weight of 1.0, this is simply the count of residential cells within that radius.

For a population density spread (smearing point-like densities over a neighbourhood):

attribute<float32> pop_density_grid (land_use) := ...;  // population per cell

// Spread population over a 1km radius using distance-decay weights
unit<spoint> kernel_1km := range(spoint, point(-10s, -10s), point(11s, 11s))
{
   attribute<uint32>  dist2_val := dist2(point(0s, 0s, kernel_1km), uint32);
   attribute<float32> weight    := dist2_val <= 100 ? 1.0f / float32(max_elem(dist2_val, 1)) : 0.0f;
}

attribute<float32> pop_potential (land_use) :=
    potential(pop_density_grid, kernel_1km/weight);

Performance: the potential function uses Fast Fourier Transformation (FFT) internally. This makes it much faster than naive nested-loop implementations, even for large kernels. A 50km kernel on a national 1km grid that would take hours in a GIS raster tool can compute in minutes.

proximity: maximum value in neighbourhood

proximity(grid_attribute, kernel_weights) is similar to potential but computes the maximum weighted value in the neighbourhood rather than the sum. This is useful for reachability questions: “is there at least one facility of type X within distance Y?”

// Is there at least one train station within 2km of each cell?
attribute<uint8> has_station_grid (land_use) := ...; // 1 where a station exists, 0 elsewhere

unit<spoint> kernel_2km := range(spoint, point(-20s, -20s), point(21s, 21s))
{
   attribute<uint32>  dist2_val := dist2(point(0s, 0s, kernel_2km), uint32);
   attribute<float32> weight    := dist2_val <= 400 ? 1.0f : 0.0f;
}

attribute<float32> station_within_2km (land_use) :=
    proximity(float32(has_station_grid), kernel_2km/weight);
// Result > 0 means at least one station is within 2km

diversity: count of distinct values in neighbourhood

diversity(grid_attribute, radius, is_circle) counts how many distinct values occur within the neighbourhood of each cell. Unlike potential and proximity, it does not take a kernel: the neighbourhood is defined by a radius in cells (a uint16 parameter) and a flag indicating a circular (1w) or square (0w) neighbourhood. This is useful for measuring land-use heterogeneity or habitat fragmentation:

// How many different land-use classes occur within 500m (5 cells) of each cell?
attribute<uint8> landuse_diversity (land_use) :=
    diversity(land_use/GridData, 5w, 1w);

The result has the same values unit as the input grid.

connected region analysis: district

district(grid_attribute) identifies connected regions: groups of adjacent cells (horizontally and vertically, not diagonally) that share the same value. The result is a new domain unit, and a Districts attribute maps each cell to its region index.

// Label connected patches of each land-use type
unit<uint32> land_use_patches := district(land_use/GridData)
{
   // Districts attribute: for each grid cell, which patch does it belong to?
   // (automatically generated as a subitem)
}

// Count cells per patch
attribute<uint32> patch_size (land_use_patches) :=
    sum(const(1[uint32], land_use), land_use_patches/Districts);

Note: null cells are treated as a special value and grouped together. To avoid merging all nodata cells into one region, filter them out before applying district.

For diagonal adjacency (8-connectivity instead of 4-connectivity), use district_8:

unit<uint32> patches_8connected := district_8(land_use/GridData);

aggregating from grid to vector

Once you have a relation from the grid to a vector domain (from poly2grid), you can aggregate grid statistics back to the vector domain using the standard aggregation functions from Module 1c:

// Rasterise COROP regions
attribute<CBS_COROP> corop_grid (land_use) :=
    poly2grid(CBS_COROP/geometry, land_use);

// Count total cells per COROP (proxy for area)
attribute<uint32> n_cells (CBS_COROP) :=
    sum(const(1[uint32], land_use), corop_grid);

// Count residential cells per COROP
attribute<uint32> n_residential (CBS_COROP) :=
    sum(uint32(land_use/GridData == 1b), corop_grid);

// Share of residential cells
attribute<float32> pct_residential (CBS_COROP) :=
    float32(n_residential) / float32(n_cells);

This pattern (rasterise, compute a grid statistic, aggregate back to vector) is one of the most common workflows in land-use and accessibility modelling.

try it yourself! [TO BE ADDED!]

  • Download the project here if you have not done so already, and unzip it to C:/prj/GeoDMSAcademy.
  • Open exercise.dms (in the GeoDMS_Academy/basic_analyses_grid/cfg subfolder) in your text editor. The configuration already includes the land-use GeoTiff and COROP shapefile from Module 2A.

Work through the following exercises:

  1. Rasterise: use poly2grid to rasterise the COROP polygons onto the land-use grid. Visualise the result in the GUI by colouring the grid cells by COROP index. How many cells fall outside all COROP regions (null values)?
  2. Aggregate: for each COROP region, compute the share of cells with land-use class 1 (residential). Which COROP has the highest residential share?
  3. Neighbourhood potential: define a circular kernel with a 500m radius (5 cells at 100m resolution). Use potential to compute the number of residential cells within 500m of each grid cell. Visualise the result as a continuous map.
  4. Diversity: compute the land-use diversity within a 300m neighbourhood (3-cell radius). Which areas show the highest diversity?
  5. District: identify connected patches of a single land-use class of your choice. Compute the area of each patch (count of cells multiplied by cell area). What is the largest contiguous patch? Check your units: open the patch area attribute in the Detail Pages and verify that its values unit carries an area metric (m² or km²). If the metric is missing, you multiplied by a bare number instead of the cell area, and your largest patch is really a cell count in disguise.

The reference solution is in result.dms in the same cfg subfolder.

Work-in-progress Future topics:

  • element-by-element operations
  • defining and using a compacted domain and basegrid to do less work
  • use combine(raster, kernel) and third_rel := first_rel + second_rel to define focal operations.
  • Grid2Polygon: https://www.geodms.nl/docs/grid-2-polygon-example.html
  • grid2polygon scripting to simplify zoning data while maintaining topological relations: https://github.com/ObjectVision/GeoDMS/wiki/Poly-to-grid-to-(simplified)-polygon

Go to previous module: Module 4, Basic analyses with vector data

Go to next module: Module 6a, Working with networks over a road network