PolarDB:Areal precipitation analysis based on GanosBase raster engine

Data model

A GanosBase raster dataset is organized as follows:

• Raster: A raster dataset — for example, a remote sensing image or a TIFF file.

• Tile: The basic storage unit of a raster object. Each tile is 256 × 256 pixels by default.

• Band: A single 2D layer in a raster dataset. Each band consists of multiple tiles, each with its own coordinate numbers.

• Cell: A pixel in a tile. Cells support multiple data types, including byte, short, int, and double.

• Pyramid: A series of reduced-resolution versions of a raster used to speed up display. Level 0 is the original data; each higher level is a coarser resolution.

• Metadata: Raster metadata including spatial extent, projection type, and pixel type.

In the database, a raster image is stored as a raster object. Each raster object is logically divided into bands, each holding the raw image data from the import process. Storage and management happen at the tile level. Tile size defaults to 256 × 256 pixels but can be changed. Each tile contains one or more bands, which in turn consist of cells representing pixels. Each raster object has its own metadata — spatial extent, data type, projection information, and coordinate numbers. If a raster is organized as a pyramid, each band has its own pyramid levels.

Advantages

GanosBase Raster uses an object-oriented storage structure: each raster — whether a remote sensing image or a DEM file — is stored as a single database record with a one-to-one mapping. Each record can hold a raster up to 1 TB. Direct operations on tiles are not permitted, which protects metadata integrity. The extension also integrates natively with time series data.

For storage efficiency, raster metadata is kept in the database while raster attribute data is stored on OSS (Object Storage Service), a low-cost object storage backend. This architecture makes large-scale raster analysis — such as processing thousands of images — practical at significantly lower storage costs.

Beyond standard spatial operations, GanosBase Raster provides statistical and algebraic operators, image color-balancing algorithms, and accelerated overview map rendering, all optimized for large-scale raster datasets.

Generate contour lines

The following statement generates contour lines at an interval of 1.0 unit, clips them to the administrative polygon, and stores the result in rs_contours:

CREATE TABLE rs_contours AS
SELECT id, max_value, min_value, ST_Intersection(a.geom,p.geom) FROM
(SELECT (ST_Contour(rast,1,'{"interval":"1.0"}')).*
FROM raster_table WHERE id=1) a, polygon_table AS p;

Import rs_contours into QGIS for an overlay display. In the image below, green numbers are measured values at observation sites and red numbers are contour values:

Generate filled surfaces

Add "polygonize":"true" to produce filled polygon surfaces instead of lines. Each output polygon covers the region between two adjacent contour values:

CREATE TABLE rs_contours_polygon AS
SELECT id, max_value, min_value, ST_Intersection(a.geom,p.geom) FROM
(SELECT (ST_Contour(rast,1,'{"interval":"1.0","polygonize":"true"}')).*
FROM raster_table WHERE id=1) a, polygon_table AS p;

The following image shows the filled surfaces colored and overlaid in QGIS:

Step 1: Review zonal statistics

Before pixel-level traversal, use ST_Statistics to get an overview of precipitation distribution within the target polygon. This function accepts a polygon geometry and a value range string, and returns per-interval statistics.

The (0,5,10,15,20] range string defines four half-open intervals: (0–5], (5–10], (10–15], and (15–20].

SELECT (ST_Statistics(rast, ST_GeomFromText('POLYGON((119.0969 28.0519, 118.9058 27.8942, 119.0502 27.5649, 119.3347 27.6292, 119.4262 27.8775, 119.4927 28.1823, 119.3812 28.1186, 119.0969 28.0519))'),
  0, '(0,5,10,15,20]',false)).*
FROM raster_table
WHERE id=1;

Expected output:

name     | band |        min         |        max         |        mean        |        sum         | count |        std         |       median       |        mode
---------+------+--------------------+--------------------+--------------------+--------------------+-------+--------------------+--------------------+--------------------
 full    |    0 | 2.5009584426879883 | 16.841625213623047 |  6.618387373747887 | 243748.58858776093 | 36829 | 2.6857099819905033 |  6.045845985412598 |  3.892089605331421
 (0-5]   |    0 | 2.5009584426879883 |  4.999907970428467 |  4.144016098134749 | 50963.109974861145 | 12298 | 0.5534420165252167 |  4.231814861297607 |  3.892089605331421
 (5-10]  |    0 |  5.000039577484131 |  9.999935150146484 |  6.852788502446667 |  136110.0852355957 | 19862 | 1.2932534814382863 |  6.623260498046875 |  5.028009414672852
 (10-15] |    0 |  10.00007438659668 | 14.993864059448242 | 11.947531793007881 |  52999.25103378296 |  4436 | 1.2279701295194279 | 11.900737762451172 | 12.430845260620117
 (15-20] |    0 |  15.00446891784668 | 16.841625213623047 | 15.777434950734413 |  3676.142343521118 |   233 | 0.5092360295014834 |   15.7352294921875 |  15.00446891784668
(5 rows)

From this output: the maximum precipitation in the area is 16.84 mm, the minimum is 2.5 mm, and there are 4,436 pixels with precipitation amounts greater than 10 mm but less than or equal to 15 mm.

Step 2: Clip the raster to the target polygon

Clip the raster to the target polygon and store the result as id=3:

INSERT INTO raster_table
SELECT 3, ST_ClipToRast(rast, ST_GeomFromText('POLYGON((119.0969 28.0519, 118.9058 27.8942, 119.0502 27.5649, 119.3347 27.6292, 119.4262 27.8775, 119.4927 28.1823, 119.3812 28.1186, 119.0969 28.0519))'),
  0, '', NULL, '', '{"chunktable":"rbt"}')
FROM raster_table WHERE id=1;

Get the pixel dimensions of the clipped raster to determine the traversal bounds:

SELECT ST_Width(rast), ST_Height(rast) FROM raster_table WHERE id=3;

Expected output:

st_width  | st_height
----------+-----------
      185 |       217

Step 3: Get the area and precipitation value for each pixel

Create a table to store per-pixel data:

CREATE TABLE pixel_pre
(
   row integer,
   col integer,
   geom geometry, -- pixel area
   pre  double precision -- pixel value
);

UPDATE raster_table SET rast=ST_SetSrid(rast, 4326);

-- Traverse each pixel of the clipped raster to get its covered area and precipitation value
DO
$do$
BEGIN
   FOR i IN 0..184 LOOP
      FOR j IN 0..216 LOOP
            INSERT INTO pixel_pre
            SELECT i,j, ST_PixelAsPolygon(rast,j,i),ST_Value(rast,0,i,j) as value FROM raster_table where id=3;
      END LOOP;
   END LOOP;
END
$do$;

This pixel-level loop is necessary because it computes the precise geographic area of each pixel polygon individually. ST_Statistics returns aggregate statistics but does not provide per-pixel area weights needed for the final calculation.

Verify the pixel data:

SELECT ROW,col,ST_AsText(geom),pre FROM pixel_pre LIMIT 10;

Expected output (each row shows the row/column index, the polygon geometry of the pixel, and its precipitation value):

row  | col |                                                                                    st_astext                                                                                     |        pre
-----+-----+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+--------------------
 107 |  33 | POLYGON((119.244752593338 28.0907410103828,119.244752593338 28.0878755431622,119.2479422912 28.0878755431622,119.2479422912 28.0907410103828,119.244752593338 28.0907410103828)) |   5.66606330871582
 107 |  34 | POLYGON((119.244752593338 28.0878755431622,119.244752593338 28.0850100759417,119.2479422912 28.0850100759417,119.2479422912 28.0878755431622,119.244752593338 28.0878755431622)) |  5.698885917663574
 107 |  35 | POLYGON((119.244752593338 28.0850100759417,119.244752593338 28.0821446087211,119.2479422912 28.0821446087211,119.2479422912 28.0850100759417,119.244752593338 28.0850100759417)) | 5.7316670417785645
 107 |  36 | POLYGON((119.244752593338 28.0821446087211,119.244752593338 28.0792791415006,119.2479422912 28.0792791415006,119.2479422912 28.0821446087211,119.244752593338 28.0821446087211)) |  5.764394283294678
 107 |  37 | POLYGON((119.244752593338 28.0792791415006,119.244752593338 28.0764136742801,119.2479422912 28.0764136742801,119.2479422912 28.0792791415006,119.244752593338 28.0792791415006)) |  5.797055244445801
 107 |  38 | POLYGON((119.244752593338 28.0764136742801,119.244752593338 28.0735482070595,119.2479422912 28.0735482070595,119.2479422912 28.0764136742801,119.244752593338 28.0764136742801)) |  5.829642295837402
 107 |  39 | POLYGON((119.244752593338 28.0735482070595,119.244752593338 28.070682739839,119.2479422912 28.070682739839,119.2479422912 28.0735482070595,119.244752593338 28.0735482070595))   |  5.862148761749268
 107 |  40 | POLYGON((119.244752593338 28.070682739839,119.244752593338 28.0678172726184,119.2479422912 28.0678172726184,119.2479422912 28.070682739839,119.244752593338 28.070682739839))    |  5.894571304321289
 107 |  41 | POLYGON((119.244752593338 28.0678172726184,119.244752593338 28.0649518053979,119.2479422912 28.0649518053979,119.2479422912 28.0678172726184,119.244752593338 28.0678172726184)) |   5.92690896987915
 107 |  42 | POLYGON((119.244752593338 28.0649518053979,119.244752593338 28.0620863381773,119.2479422912 28.0620863381773,119.2479422912 28.0649518053979,119.244752593338 28.0649518053979)) |  5.959161758422852
(10 rows)

Step 4: Calculate weighted areal precipitation

Raster pixels in EPSG:4326 (geographic coordinates) have areas that vary with latitude, so convert each pixel polygon to EPSG:3857 (metric projection) before computing its area in square meters. Then compute the area-weighted average:

SELECT sum(ST_Area(ST_Transform(geom,3857)) * pre) / sum(ST_Area(ST_Transform(geom,3857))) as "precipitation(mm/m^2)"
FROM pixel_pre WHERE pre>0;

Expected output:

precipitation(mm/m^2)
-----------------------
     9.434021577002174
(1 row)

The areal precipitation for this polygon is approximately 9.43 mm/m².