Créer une table raster carrelée vide dans PostGIS

Question

Je suis nouveau à postgres/postgis et j'essaie de comprendre comment faire un raster carrelé vide avec postgis. Je voudrais générer un raster vide contenant 5000 x 2000 cellules, que je voudrais interroger plus tard pour trouver une cellule spécifique à l'emplacement x/y ou ajouter une valeur de cellule à l'emplacement x/y. J'ajouterai fondamentalement des valeurs de cellules individuelles au raster vide à l'avenir (par exemple, signaler des observations d'animaux dans des cellules au fil du temps).

Je trouve que la plupart de ce que cela est possible en créant d'abord une table raster:

CREATE TABLE myRasterTable(rid serial primary key, rast raster); 

, puis en ajoutant une trame vide:

INSERT INTO myRasterTable(rid,rast) 
VALUES(1, ST_MakeEmptyRaster(5000, 2000, 2485869.5728, 1299941.7864, 100, 100, 0, 0, 2056)); 

(Référence: http://suite.opengeo.org/docs/latest/dataadmin/pgGettingStarted/raster2pgsql.html)

J'ai également ajouté un groupe et ajouté une valeur à l'une des cellules raster:

UPDATE myRasterTable SET rast = ST_AddBand(rast, 1, '32BF'::text, 0); 
UPDATE myRasterTable SET rast = ST_SetValue(rast, 1,ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056),987.654321) 

(Référence: https://gis.stackexchange.com/questions/14960/postgis-raster-value-of-a-lat-lon-point)

Je peux maintenant interroger la valeur de la cellule raster j'ai mis:

// Location with value 
SELECT rid, ST_Value(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056),false) val FROM myRasterTable 
// Return = 987.654296875 in 0.5224609375 Seconds 

// Location without value: 
SELECT rid, ST_Value(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.0,48.5),4326),2056),false) val FROM myRasterTable 
// Return = 0 in 0.51311993598938 Seconds 

J'ai trouvé de nombreux messages disant que pour les rasters plus importants, il est essentiel ce qui a trait à la performance de requête, que rasters sont carrelées avec des degrés d'environ 100X100 cellules par tuile (références: http://postgis.17.x6.nabble.com/raster-loading-and-ST-Value-performance-td4999924.html ET https://duncanjg.wordpress.com/2013/09/21/effect-of-tile-size-and-data-storage-on-postgis-raster-query-times/)

maintenant, je ne sais pas:

1. Comment créer des carreaux à partir de mon grand raster dans PostGIS ou comment créer un carreau matriciel dès le départ avec PostGIS? En d'autres termes: Quelle requête dois-je utiliser pour créer une collection de carreaux raster couvrant un raster 5000X2000 avec des carreaux d'une étendue de 100X100 cellules contenant plusieurs bandes?
2. Une fois les tuiles créées, ai-je besoin de créer un index spatial séparé ou est-ce fait automatiquement?
3. Enfin, comment puis-je effectuer une requête de quelle valeur de cellule a un emplacement spécifique après que le raster a été mis en mosaïque?

Toute aide serait appréciée!

Answer 1

Après un certain nombre d'essais et d'erreurs, il semble que j'ai trouvé la solution.

Pour créer une nouvelle table dans PostGIS contenant des tuiles vides qui ajoutent une extension de 5000 X 2000 cellules, je cette requête:

CREATE TABLE myRasterTable (rid integer, rast raster); 
INSERT INTO myRasterTable(rast) VALUES(ST_Tile(ST_MakeEmptyRaster(5000, 2000, 2485869.5728, 1299941.7864, 100, 100, 0, 0, 2056), 10,10, TRUE, NULL)); 

Cela fait une trame temporaire avec une étendue de 5000 X 2000 en utilisant ST_MakeEmptyRaster. Il colmate ensuite le raster temporaire en utilisant ST_Tile dans des carreaux de cellule 10X10 et ajoute ces carreaux à la table.

Je peux alors ajouter mon groupe:

UPDATE myRasterTable SET rast = ST_AddBand(rast, 1, '32BUI'::text, 0, NULL); 

Et enfin, je peux ajouter et récupérer les valeurs avec:

// Add cell value 
UPDATE myRasterTable SET rast = ST_SetValue(rast, 1,ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056),987654321) WHERE ST_Intersects(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056)); 

// Get cell value 
SELECT rid, ST_Value(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056), false) FROM myRasterTable WHERE ST_Intersects(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056)); 
// Return = 987654321 in 0.22717499732971 Seconds 

Comme vous pouvez le voir, ce déjà réduit le temps de requête à ~ 2/5 de l'heure d'origine.

Si j'ajoute en outre un indice:

CREATE INDEX myRasterTable_rast_gist_idx ON myRasterTable USING GIST (ST_ConvexHull(rast)); 

Et puis recommencez la requête, je reçois:

// Get cell value 
SELECT rid, ST_Value(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056), false) FROM myRasterTable WHERE ST_Intersects(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056)); 
// Return = 987654321 in 0.084153890609741 Seconds 

Comme vous pouvez le voir, le querytime a de nouveau coupé un peu plus de la moitié, résultant en un querytime de moins de 1/5 de la requête originale. J'ai également essayé différentes tailles de carreaux et constaté qu'une taille de carreau de 10 X 10 donnait les meilleures performances, mais seulement légèrement mieux qu'une taille de carreaux de 100 X 100.

Si quelqu'un a d'autres optimisations, n'hésitez pas à ajouter .

EDIT (08.07.2016)

J'ai écrit un petit billet de blog à ce sujet. Si vous êtes intéressé, vous pouvez jeter un oeil ici: http://www.geonet.ch/postgres-postgis-of-rasters-and-geojsons/