30. Raster¶
PostGIS stöder en annan typ av spatial datatyp som kallas raster. Rasterdata använder, precis som geometridata, kartesiska koordinater och ett spatialt referenssystem. Men istället för vektordata representeras rasterdata som en n-dimensionell matris som består av bildpunkter och band. Banden definierar antalet matriser du har. Varje pixel lagrar ett värde som motsvarar varje band. Så en 3-bandsraster, t.ex. en RGB-bild, skulle ha 3 värden för varje pixel som motsvarar banden rött-grönt-blått.
Även om vackra bilder, som de du ser på din TV-skärm, är raster, är raster kanske inte så spännande att titta på. I ett nötskal är ett raster en matris, fäst i ett koordinatsystem, som har värden som kan representera vad du vill att de ska representera.
Eftersom raster lever i kartesiskt rum kan raster interagera med geometrier. PostGIS erbjuder många funktioner som tar både raster och geometrier som indata. Många operationer som tillämpas på raster kommer att resultera i geometrier. Vanliga sådana är ST_Polygon, ST_Envelope, ST_ConvevexHull och ST_MinConvexHull enligt nedan. När du kastar ett raster till en geometri, är det som matas ut ST_ConvevexHull av rastret.
Rasterformatet används ofta för att lagra höjddata, temperaturdata, satellitdata och tematiska data som representerar saker som miljöföroreningar, befolkningstäthet och förekomster av miljöfaror. Du kan använda raster för att lagra alla numeriska data som har en meningsfull koordinatplacering. Den enda begränsningen är att de numeriska datatyperna måste vara desamma för alla data i ett visst band.
Även om rasterdata kan skapas från grunden i PostGIS är det vanligare att ladda rasterdata från olika format med hjälp av kommandot raster2pgsql som medföljer PostGIS. Innan allt detta måste du aktivera rasterstöd i din databas genom att köra kommandot:
CREATE EXTENSION postgis_raster;
30.1. Skapa Raster från Geometrier¶
Vi börjar med att skapa rasterdata från vektordata och går sedan vidare till den mer spännande metoden att ladda data från en rasterkälla. Du kommer att upptäcka att rasterdata finns tillgängliga i överflöd och ofta gratis från olika myndighetssidor.
Vi börjar med att konvertera några geometrier till raster med hjälp av funktionen ST_AsRaster på följande sätt.
CREATE TABLE rasters (name varchar, rast raster);
INSERT INTO rasters(name, rast)
SELECT f.word, ST_AsRaster(geom, width=>150, height=>150)
FROM (VALUES ('Hello'), ('Raster') ) AS f(word)
, ST_Letters(word) AS geom;
CREATE INDEX ix_rasters_rast
ON rasters USING gist(ST_ConvexHull(rast));
I exemplet ovan CREATEs en tabell (rasters) från geometrier som bildats av bokstäver med hjälp av PostGIS 3.2+ ST_Letters-funktionen. Raster kan på samma sätt som geometrier dra nytta av spatiala index. Det spatiala indexet som används för raster är ett funktionellt index som indexerar rastrets geometriska konvexskal.
Du kan se några användbara metadata för dina raster med följande fråga som använder postgis raster ST_Count-funktionen för att räkna antalet bildpunkter som har data och ST_MetaData-funktionen för att tillhandahålla alla typer av användbar bakgrundsinformation för våra raster.
SELECT name, ST_Count(rast) As num_pixels, md.*
FROM rasters, ST_MetaData(rast) AS md;
name | num_pixels | upperleftx | upperlefty | width | height | scalex | scaley | skewx | skewy | srid | numbands
--------+------------+------------+-------------------+-------+--------+--------------------+---------------------+-------+-------+------+----------
Hello | 13926 | 0 | 77.10000000000001 | 150 | 150 | 1.226888888888889 | -0.5173333333333334 | 0 | 0 | 0 | 1
Raster | 11967 | 0 | 75.4 | 150 | 150 | 1.7226319023207244 | -0.5086666666666667 | 0 | 0 | 0 | 1
(2 rows)
Observera
Det finns två nivåer av rasterfunktioner. Det finns funktioner som ST_MetaData som arbetar på rasternivå och det finns funktioner som ST_Count-funktionen och ST_BandMetaData-funktionen som arbetar på bandnivå. De flesta funktioner i postgis raster som arbetar på bandnivå arbetar bara med ett band i taget och antar att det band du vill ha är 1.
Om du har ett raster med flera band och du behöver räkna de bildpunkter som inte har några datavärden i ett annat band än 1, måste du uttryckligen ange bandnumret på följande sätt ST_Count(rast,2).
Observera att alla raster har en dimension på 150x150. Detta är inte idealiskt. Det innebär att våra raster, för att tvinga fram det, är klämda på alla möjliga sätt. Om vi bara kunde se hur fula rastren är framför oss.
30.2. Läs in raster med hjälp av raster2pgsql¶
raster2pgsql är ett kommandoradsverktyg som ofta paketeras med PostGIS. Om du är på Windows och använde applikationsstackbuilder PostGIS Bundle, hittar du raster2pgsql.exe i mappen C: \ Program Files \ PostgreSQL \ 15 \ bin där 15 ska ersättas med den version av PostgreSQL du kör.
Om du använder Postgres.app hittar du raster2pgsql bland de andra Postgres.app CLI Tools.
På Ubuntu och Debian behöver du
apt install postgis
för att ha PostGIS kommandoradsverktyg installerade. Detta kan också installera en ytterligare version av PostgreSQL. Du kan se en lista över kluster i Debian/Ubuntu med hjälp av kommandot pg_lsclusters och ta bort dem med kommandot pg_dropcluster.
För denna och senare övningar kommer vi att använda nyc_dem.tif som finns i filen PG Raster Workshop Dataset https://postgis.net/stuff/workshop-data/postgis_raster_workshop.zip. För vissa geometri/rasterexempel kommer vi också att använda NYC-data som laddats från tidigare kapitel. I stället för att ladda tif-filen kan du återställa nyc_dem.backup som ingår i zip-filen i din databas med hjälp av kommandoradsverktyget pg_restore eller menyn pgAdmin Restore.
Observera
Denna rasterdata hämtades från NYC DEM 1-foot Integer som är en 3GB DEM tif som representerar höjd i förhållande till havsnivån med byggnader och övervatten borttagna. Vi skapade sedan en lägre res-version av den.
Verktyget raster2pgsql liknar shp2gpsql förutom att istället för att ladda ESRI shapefiler till PostGIS geometri/geografitabeller, laddar det alla GDAL-stödda rasterformat till rastertabeller. Precis som i shp2pgsql kan du skicka ett SRID (spatial reference id) för källan. Till skillnad från shp2pgsql kan den härleda källdatans spatiala referenssystem om källdata har lämpliga metadata.
För en fullständig beskrivning av alla möjliga alternativ, se raster2pgsql options.
Några andra anmärkningsvärda alternativ som raster2pgsql erbjuder som vi inte kommer att täcka är:
Möjlighet att ange källans SRID. Istället förlitar vi oss på raster2pgsql gissningsförmåga.
Möjlighet att ställa in värdet nodata, om det inte anges försöker raster2pgsql härleda från filen.
Möjlighet att ladda raster utanför databasen.
För att ladda alla tif-filer i vår mapp och även skapa översikter, skulle vi köra nedanstående.
raster2pgsql -d -e -l 2,3 -I -C -M -F -Y -t 256x256 *.tif nyc_dem | psql -d nyc
-d för att ta bort tabellerna om de redan finns
Ovanstående kommando använder -e för att ladda omedelbart istället för att binda i en transaktion
-C ställer in rasterbegränsningar, detta är användbart för raster_columns för att visa information. Du kanske vill kombinera med -x för att utesluta extent-begränsningen, som är en långsam begränsning att kontrollera och som också hindrar framtida belastningar i tabellen.
-M för att dammsuga och analysera efter laddning, för att förbättra statistiken i frågeplaneraren
-Y för att använda kopiering i omgångar om 50. Om du kör PostGIS 3.3 eller senare kan du använda -Y 1000 för att kopiera i satser om 1000, eller ännu högre antal. Detta går snabbare, men använder mer minne.
-l 2,3 för att skapa överskådliga tabeller: o_2_ncy_dem och o_3_nyc_dem. Detta är användbart för visning av data.
-I för att skapa ett spatialt index
-F för att lägga till filnamn, om du bara har en tif-fil är det här ganska meningslöst. Om du hade flera, detta skulle vara användbart för att berätta vilken fil varje rad kom från.
-t för att ställa in blockstorleken. Observera att om du inte är säker på vilken storlek som är bäst att använda, använd -t auto istället och raster2pgql kommer att använda samma tiling som vad som fanns i tif. Utmatningen kommer att berätta vilken blockstorlek den valde. Avbryt om det ser stort eller konstigt ut. Originalfilen hade en storlek på 300x7 vilket inte är idealiskt.
Använd psql för att köra den genererade sql:en mot databasen. Om du vill dumpa till en fil istället, använd > nyc_dem.sql
I det här exemplet har vi bara en tif-fil, så vi kan i stället ange det fullständiga filnamnet i stället för *.tif. Om filerna inte finns i din aktuella katalog kan du också ange en mappsökväg med *.tif.
Observera
Om du använder Windows och behöver hänvisa till mappen, se till att inkludera enhetsbeteckningen, t.ex. C:/workshop/*.tif
I PostGIS-lingo används ofta termerna raster tile och raster omväxlande. En rasterplatta motsvarar egentligen ett visst raster i en rasterkolumn som är en delmängd av ett större raster, till exempel NYC dem-data som vi just laddade. Detta beror på att när raster laddas in i PostGIS från stora rasterfiler, hackas de i många rader för att göra dem hanterbara. Varje raster i varje rad är då en del av ett större raster. Varje kakel täcker samma storlek på området som anges av den blockstorlek du angav. Rasters är tyvärr begränsade av 1 GB PostgreSQL TOAST-gränsen och även den långsamma processen med detoasting och så måste vi hugga upp för att uppnå anständig prestanda eller till och med lagra dem.
30.3. Visa raster i webbläsaren¶
Även om pgAdmin och psql ännu inte har någon mekanism för att visa postgis-raster, har vi ett par alternativ. För små raster är det enklaste att mata ut till ett webbvänligt rasterformat som PNG med hjälp av batterier som ingår i postgis rasterfunktioner som ST_AsPNG eller ST_AsGDALRaster som listas i PostGIS Rasterutmatningsfunktioner. När dina raster blir större kommer du att vilja gå över till ett verktyg som QGIS för att visa dem i all sin prakt eller GDAL-familjen med kommandoradsverktyg som gdal_translate för att exportera dem till andra rasterformat. Kom dock ihåg att Postgis raster är byggda för analys, inte för att generera vackra bilder som du kan titta på.
En varning, som standard är alla olika utmatningar för rastertyper inaktiverade. För att kunna använda dessa måste du aktivera drivrutiner, alla eller en delmängd enligt beskrivningen i Aktivera GDAL Raster-drivrutiner
SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
Om du inte vill behöva göra detta för varje anslutning kan du ställa in på databasnivå med hjälp av:
ALTER DATABASE nyc SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
Varje ny anslutning till databasen kommer att använda den inställningen.
Kör nedanstående fråga och kopiera och klistra in resultatet i adressfältet i din webbläsare.
SELECT 'data:image/png;base64,' ||
encode(ST_AsPNG(rast),'base64')
FROM rasters
WHERE name = 'Hello';
För de raster som skapats hittills har vi inte specificerat antalet band eller ens definierat deras relation till jorden. Därför har våra raster ett okänt spatialt referenssystem (0).
Man kan tänka sig ett rasters exoskelett som en geometri. En matris innesluten i ett geometriskt kuvert. För att kunna göra användbara analyser måste vi georeferera våra raster, vilket innebär att vi vill att varje pixel (rektangel) ska representera någon meningsfull del av rymden.
ST_AsRaster har många överladdade representationer. Det tidigare exemplet använde den enklaste implementeringen och accepterade standardargumenten som är 8BUI och 1 band, där inga data är 0. Om du behöver använda de andra varianterna bör du använda anropssyntaxen för namngivna argument så att du inte av misstag hamnar i fel variant av funktionen eller får funktionen är inte unik-fel.
Om du börjar med en geometri som har ett spatialt referenssystem kommer du att sluta med ett raster med samma spatiala referenssystem. I nästa exempel placerar vi våra ord i New York i ljusa, glada färger. Vi kommer också att använda pixelskala i stället för bredd och höjd så att våra rasterpixelstorlekar representerar 1 meter x 1 meter utrymme.
INSERT INTO rasters(name, rast)
SELECT f.word || ' in New York' ,
ST_AsRaster(geom,
scalex => 1.0, scaley => -1.0,
pixeltype => ARRAY['8BUI', '8BUI', '8BUI'],
value => CASE WHEN word = 'Hello' THEN
ARRAY[10,10,100] ELSE ARRAY[10,100,10] END,
nodataval => ARRAY[0,0,0], gridx => NULL, gridy => NULL
) AS rast
FROM (
VALUES ('Hello'), ('Raster') ) AS f(word)
, ST_SetSRID(
ST_Translate(ST_Letters(word),586467,4504725), 26918
) AS geom;
Om vi sedan tittar på detta ser vi en geometri som inte är färgad med kvadrater.
SELECT 'data:image/png;base64,' ||
encode(ST_AsPNG(rast),'base64')
FROM rasters
WHERE name = 'Hello in New York';
Upprepa för Raster:
SELECT 'data:image/png;base64,' ||
encode(ST_AsPNG(rast),'base64')
FROM rasters
WHERE name = 'Raster in New York';
Vad som är mer talande är att om vi kör om
SELECT name, ST_Count(rast) As num_pixels, md.*
FROM rasters, ST_MetaData(rast) AS md;
Observera metadata för New York-posterna. De har det spatiala referenssystemet New York state plane meter. De har också samma skala. Eftersom varje enhet är 1x1 meter är bredden på ordet Raster nu bredare än Hello.
name | num_pixels | upperleftx | upperlefty | width | height | scalex | scaley | skewx | skewy | srid | numbands
-------------------+------------+------------+-------------------+-------+--------+--------------------+---------------------+-------+-------+-------+----------
Hello | 13926 | 0 | 77.10000000000001 | 150 | 150 | 1.226888888888889 | -0.5173333333333334 | 0 | 0 | 0 | 1
Raster | 11967 | 0 | 75.4 | 150 | 150 | 1.7226319023207244 | -0.5086666666666667 | 0 | 0 | 0 | 1
Hello in New York | 8786 | 586467 | 4504802.1 | 184 | 78 | 1 | -1 | 0 | 0 | 26918 | 3
Raster in New York | 10544 | 586467 | 4504800.4 | 258 | 76 | 1 | -1 | 0 | 0 | 26918 | 3
(4 rows)
30.4. Raster Spatiala Katalogtabeller¶
I likhet med geometri- och geograftyperna har raster en uppsättning kataloger som visar alla rasterkolumner i din databas. Dessa är raster_columns och raster_overviews.
30.4.1. raster_columns¶
Vyn raster_columns är syskon till geometry_columns och geography_columns och tillhandahåller i stort sett samma data och mer därtill, men för rasterkolumner.
SELECT *
FROM raster_columns;
Om du utforskar tabellen hittar du det här:
r_table_catalog | r_table_schema | r_table_name | r_raster_column | srid | scale_x | scale_y | blocksize_x | blocksize_y | same_alignment | regular_blocking | num_bands | pixel_types | nodata_values | out_db | extent | spatial_index
----------------+----------------+--------------+-----------------+------+---------+---------+-------------+-------------+----------------+------------------+-----------+-------------+---------------+--------+--------+---------------
nyc | public | rasters | rast | 0 | | | | | f | f | | | | | | t
nyc | public | nyc_dem | rast | 2263 | 10 | -10 | 256 | 256 | t | f | 1 | {16BUI} | {NULL} | {f} | | t
nyc | public | o_2_nyc_dem | rast | 2263 | 20 | -20 | 256 | 256 | t | f | 1 | {16BUI} | {NULL} | {f} | | t
nyc | public | o_3_nyc_dem | rast | 2263 | 30 | -30 | 256 | 256 | t | f | 1 | {16BUI} | {NULL} | {f} | | t
(4 rows)
en nedslående rad med i stort sett ofylld information för rasters tabell.
Till skillnad från geometri och geografi har raster inte stöd för typmodifierare, eftersom utrymmet för typmodifierare är för begränsat och det finns fler kritiska egenskaper än vad som ryms i en typmodifierare.
Raster förlitar sig istället på begränsningar och läser tillbaka dessa begränsningar som en del av vyn.
Titta på de andra raderna från de tabeller vi laddade med raster2pgsql. Eftersom vi använde -C-omkopplaren lade raster2pgsql till begränsningar för srid och annan information som den kunde läsa från tif eller som vi skickade in. Översiktstabellerna som genererades med -l switch o_2_nyc_dem och o_3_nyc_dem dyker också upp.
Låt oss försöka lägga till några begränsningar i vår tabell.
SELECT AddRasterConstraints('public'::name, 'rasters'::name, 'rast'::name);
Och du kommer att bombarderas med en hel massa meddelanden om hur dina rasterdata är en röra och ingenting kan begränsas. Om du tittar på raster_columns igen, är det fortfarande samma nedslående historia med många tomma rader för rasters.
För att begränsningar ska kunna tillämpas måste alla raster i tabellen vara möjliga att begränsa med minst en regel.
Vi kan kanske göra det här, låt oss bara ljuga och säga att alla våra data finns i New York State-planet.
UPDATE rasters SET rast = ST_SetSRID(rast,26918)
WHERE ST_SRID(rast) <> 26918;
SELECT AddRasterConstraints('public'::name, 'rasters'::name, 'rast'::name);
SELECT r_table_name AS t, r_raster_column AS c, srid,
blocksize_x AS bx, blocksize_y AS by, scale_x AS sx, scale_y AS sy,
ST_AsText(extent) AS e
FROM raster_columns
WHERE r_table_name = 'rasters';
Ah framsteg:
t | c | srid | bx | by | sx | sy | e
----------+------+-------+-----+-----+----+----+------------------------------------------
rasters | rast | 26918 | 150 | 150 | | | POLYGON((0 -0.90000000000..
(1 row)
Ju mer du kan begränsa alla dina raster, desto fler kolumner kommer du att se ifyllda och desto fler operationer kommer du att kunna göra över alla plattor i ditt raster. Tänk på att du i vissa fall kanske inte vill tillämpa alla begränsningar.
Om du till exempel planerar att ladda mer data i din rastertabell vill du hoppa över extent-begränsningen eftersom det skulle kräva att alla raster ligger inom extent-begränsningens omfattning.
30.4.2. raster_overviews¶
Rasteröversiktskolumner visas både i metakatalogen raster_columns och i en annan metakatalog som heter raster_overviews. Översikter används främst för att snabba upp visningen vid högre zoomnivåer. De kan också användas för snabb analys av baksidan av kuvertet, vilket ger mindre exakt statistik, men med en mycket snabbare hastighet än att tillämpa på den råa rastertabellen.
Kör för att inspektera översikterna:
SELECT *
FROM raster_overviews;
och du kommer att se resultatet:
o_table_catalog | o_table_schema | o_table_name | o_raster_column | r_table_catalog | r_table_schema | r_table_name | r_raster_column | overview_factor
----------------+----------------+--------------+-----------------+-----------------+----------------+--------------+-----------------+-----------------
nyc | public | o_2_nyc_dem | rast | nyc | public | nyc_dem | rast | 2
nyc | public | o_3_nyc_dem | rast | nyc | public | nyc_dem | rast | 3
(2 rows)
Tabellen raster_overviews ger dig bara overview_factor och namnet på den överordnade tabellen. All denna information är något som du själv kunde ha räknat ut med hjälp av raster2pgsql namnkonvention för översikter.
Overview-faktorn anger vilken upplösning raden har i förhållande till sin överordnade enhet. En overview_factor på 2 innebär att 2x2 = 4 brickor får plats på en overview_2-bricka. På samma sätt innebär en overview_factor på 3 att 2x2x2 = 8 brickor av originalet kan tryckas in i en overview_3-bricka.
30.5. Vanliga rasterfunktioner¶
Tillägget postgis_raster har över 100 funktioner att välja mellan. PostGIS rasterfunktionalitet har utformats efter PostGIS geometristöd. Du kommer att hitta en överlappning av funktioner mellan raster och geometri där det är vettigt. Vanliga funktioner som du kommer att använda och som har motsvarigheter i geometrivärlden är ST_Intersects, ST_SetSRID, ST_SRID, ST_Union, ST_Intersection och ST_Transform.
Förutom dessa överlappande funktioner har programmet stöd för överlappningsoperatorn && mellan raster och mellan ett raster och geometri. Det finns också många funktioner som fungerar tillsammans med geometri eller som är mycket specifika för raster.
Du behöver en funktion som ST_Union för att återskapa en region. Eftersom prestanda blir långsammare ju fler bildpunkter en funktion behöver analysera, behöver du en snabbverkande funktion ST_Clip för att klippa rastren till bara de delar som är intressanta för din analys.
Slutligen behöver du ST_Intersects eller && för att zooma in på de rasterplattor som innehåller dina intresseområden. Operatorn && är en snabbare process än ST_Intersects. Båda kan dra nytta av spatiala index för raster. Vi går först igenom dessa grundläggande funktioner innan vi går vidare till andra avsnitt där vi använder dem tillsammans med andra raster- och geometrifunktioner.
30.5.1. Sammanslagning av raster med ST_Union¶
Funktionen ST_Union för raster, precis som geometrins motsvarighet ST_Union, sammanför en uppsättning raster till ett enda raster. Precis som med geometri kan dock inte alla raster kombineras, men reglerna för rasterunionering är mer komplicerade än geometrireglerna. När det gäller geometrier räcker det med att ha samma spatiala referenssystem, men för raster är det inte tillräckligt.
Om du skulle försöka, följande:
SELECT ST_Union(rast)
FROM rasters;
Du skulle bli summariskt straffad med ett fel:
ERROR: rt_raster_from_two_rasters: De två raster som tillhandahålls har inte samma inriktning SQL-status: XX000
Vad är det för något med denna inriktning som hindrar dig från att sammanfoga dina värdefulla raster?
För att raster ska kunna kombineras måste de så att säga vara på samma rutnät. Det innebär att de måste ha samma pixelstorlekar, samma orientering (skevheten), samma spatiala referenssystem och att deras bildpunkter inte får korsa in i varandra, vilket innebär att de delar samma världsliga pixelnät.
Om du försöker samma fråga, men bara med ord som vi noggrant placerat i New York.
Återigen, samma fel. Det är samma spatiala ref-system, samma pixelstorlekar, och ändå är det inte tillräckligt bra. För att deras rutnät är fel.
Vi kan åtgärda detta genom att flytta de övre vänstra y-koordinaterna en aning och sedan försöka igen. Om våra rutnät börjar på heltalsnivå eftersom våra pixelstorlekar är hela heltal, kommer bildpunkterna inte att korsa in i varandra.
UPDATE rasters SET rast = ST_SetUpperLeft(rast,
ST_UpperLeftX(rast)::integer,
ST_UpperLeftY(rast)::integer)
WHERE name LIKE '%New York';
SELECT ST_Union(rast ORDER BY name)
FROM rasters
WHERE name LIKE '%New York%';
Voila det fungerade, och om vi skulle visa, skulle vi se något liknande detta:
Observera
Om du någonsin är osäker på varför dina raster inte har samma inriktning kan du använda funktionen ST_SameAlignment, som jämför 2 raster eller en uppsättning raster och talar om för dig om de har samma inriktning. Om du har aktiverat notiser kommer NOTICE att berätta vad som är fel med rastrerna i fråga. Med ST_NotSameAlignmentReason, istället för bara ett meddelande, kommer orsaken att matas ut. Det fungerar dock bara med två raster åt gången.
Ett viktigt sätt på vilket rasterfunktionen ST_Union(raster) avviker från geometrifunktionen ST_Union(geometry) är att den tillåter ett argument som kallas uniontype. Detta argument är som standard inställt på LAST om du inte anger det, vilket innebär att du tar de LAST rasterpixelvärdena vid tillfällen där rasterpixelvärdena överlappar varandra. Som en allmän regel ignoreras bildpunkter i ett band som är markerade som no-data.
Precis som med de flesta aggregat i PostgreSQL kan du lägga till en :command: ORDER BY -klausul som en del av funktionsanropet som görs i det tidigare exemplet. Genom att ange ordningen kan du styra vilken raster som prioriteras. Så i vårt tidigare exempel trumfade Raster Hello eftersom Raster är alfabetiskt sist.
Observera, om du byter ordning:
SELECT ST_Union(rast ORDER BY name DESC)
FROM rasters
WHERE name LIKE '%New York%';
Sedan övertrumfar Hello Raster eftersom Hello nu är den sista överlagrade.
Unionstypen FIRST är motsatsen till LAST.
Men ibland kanske LAST inte är rätt operation. Låt oss anta att våra raster representerade två olika uppsättningar observationer från två olika enheter. Dessa enheter mäter samma sak, och vi är inte säkra på vilken som är rätt när de korsar vägar, så vi skulle istället vilja ta ”medelvärdet” av resultaten. Vi skulle göra så här:
SELECT ST_Union(rast, 'MEAN')
FROM rasters
WHERE name LIKE '%New York%';
Voila det fungerade, och om vi skulle visa, skulle vi se något liknande detta:
Så istället för att trumfa, har vi en blandning av de två krafterna. När det gäller unionstypen MEAN är det meningslöst att ange ordning, eftersom resultatet blir ett genomsnitt av överlappande pixelvärden.
Observera att för geometrier, eftersom geometrier är vektorer och därmed inte har några värden förutom där eller inte där, finns det egentligen ingen tvetydighet om hur man kombinerar två vektorer när de korsar varandra.
En annan egenskap hos raster ST_Union som vi glömde bort är tanken på om du ska returnera alla band eller bara vissa band. När du inte anger vilka band som ska förenas kommer ST_Union att kombinera samma bandnummer och använda LAST-strategin för förening. Om du har flera band är detta kanske inte vad du vill göra. Du kanske bara vill förena det andra bandet. I det här fallet är det det gröna bandet och du vill ha antalet pixelvärden.
SELECT ST_BandPixelType(ST_Union(rast, 2, 'COUNT'))
FROM rasters
WHERE name LIKE '%New York%';
st_bandpixeltype
------------------
32BUI
(1 row)
Observera att när det gäller unionstypen COUNT, som räknar antalet bildpunkter som fyllts i och returnerar det värdet, är resultatet alltid ett 32BUI på samma sätt som när du gör ett COUNT i sql, är resultatet alltid ett bigint, för att tillgodose stora räkningar.
I andra fall ändras inte bandpixeltypen och den ställs in på maxvärdet eller avrundas om beloppen överskrider gränserna för typen. Varför skulle någon någonsin vilja räkna bildpunkter som korsar varandra på en plats. Anta att alla dina raster representerar polisskvadroner och incidenter med arresteringar i områdena. Varje värde kan representera en annan typ av arresteringsorsak. Du gör statistik över hur många arresteringar i varje region, därför bryr du dig bara om antalet arresteringar.
Eller kanske vill du göra alla band, men du vill ha olika strategier.
SELECT ST_Union(rast, ARRAY[(1, 'MAX'),
(2, 'MEAN'),
(3, 'RANGE')]::unionarg[])
FROM rasters
WHERE name LIKE '%New York%';
Om du använder varianten unionarg[] av funktionen ST_Union kan du också blanda om ordningen på banden.
30.5.2. Klippning av Rasters med hjälp av ST_Intersects¶
Funktionen ST_Clip är en av de mest använda funktionerna för PostGIS-raster. Huvudskälet är att ju fler bildpunkter du behöver inspektera eller göra operationer på, desto långsammare blir din bearbetning. ST_Clip klipper ditt raster till bara det intressanta området, så att du kan isolera dina operationer till bara det området.
Den här funktionen är också speciell eftersom den utnyttjar geometrins möjligheter för att underlätta rasteranalys. För att minska antalet bildpunkter som ST_Union måste hantera, klipps varje raster först till det område som vi är intresserade av.
SELECT ST_Union( ST_Clip(r.rast, g.geom) )
FROM rasters AS r
INNER JOIN
ST_Buffer(ST_Point(586598, 4504816, 26918), 100 ) AS g(geom)
ON ST_Intersects(r.rast, g.geom)
WHERE r.name LIKE '%New York%';
Detta exempel visar flera funktioner som arbetar tillsammans. Funktionen ST_Intersects som används är den som medföljer postgis_raster och kan korsa 2 raster eller ett raster och en geometri. I likhet med geometri ST_Intersects kan raster ST_Intersects dra nytta av spatiala index på raster- eller geometritabellerna.
Observera
Som standard kommer ST_Clip att utelämna bildpunkter där pixelns centroid inte skär geometrin. Detta kan vara irriterande för stora bildpunkter och du kanske föredrar att istället inkludera en pixel om någon del av pixeln vidrör geometrin. Argumentet touched introducerades i PostGIS 3.5. Ersätt din ST_Clip(r.rast, g.geom) med ST_Clip(r.rast, g.geom, touched => true) och voila alla bildpunkter som skär din geometri på något sätt kommer att inkluderas.
30.6. Konvertera Raster till Geometrier¶
Raster kan lika gärna omvandlas till geometrier.
30.6.1. Polygonen i ett raster med ST_Polygon¶
Låt oss börja med vårt tidigare exempel, men konvertera det till en polygon med hjälp av funktionen ST_Polygon.
SELECT ST_Polygon(ST_Union( ST_Clip(r.rast, g.geom) ))
FROM rasters AS r
INNER JOIN
ST_Buffer(ST_Point(586598, 4504816, 26918), 100 ) AS g(geom)
ON ST_Intersects(r.rast, g.geom)
WHERE r.name LIKE '%New York%';
Om du klickar på geometrivisaren i pgAdmin kan du se detta i all sin prakt utan några hack.
ST_Polygon tar hänsyn till alla bildpunkter som har värden (inte inga data) i ett visst band och omvandlar dem till geometri. Liksom många andra funktioner i raster beaktar ST_Polygon endast 1 band. Om inget band anges beaktas endast det första bandet.
30.6.2. Pixelrektanglarna i ett raster med ST_PixelAsPolygons¶
En annan ofta använd funktion är funktionen ST_PixelAsPolygons. Du bör sällan använda ST_PixelAsPolygons på ett stort raster utan att först klippa eftersom du kommer att få miljontals rader, en för varje pixel.
ST_PixelAsPolygons returnerar en tabell som består av geom, val, x och y. Där x är kolumnnumret och y är radnumret i rastret.
ST_PixelAsPolygons i likhet med andra rasterfunktioner arbetar med ett band i taget och arbetar med band 1 om inget band anges. Den returnerar också som standard endast bildpunkter som har värden.
SELECT gv.*
FROM rasters AS r
CROSS JOIN LATERAL ST_PixelAsPolygons(rast) AS gv
WHERE r.name LIKE '%New York%'
LIMIT 10;
Vilket matar ut:
och om vi inspekterar med geometrivisaren, ser vi:
Om vi vill ha alla bildpunkter i alla våra band måste vi göra något liknande nedan. Notera skillnaderna i detta exempel från tidigare.
1. Setting exclude_nodata_value to make sure all pixels are returned so that our sets of calls return the same number of rows. The rows out of the function will be naturally in the same order.
2. Using the PostgreSQL ROWS FROM constructor , and aliasing each set of columns from our function output with names. So for example the band 1 columns (geom, val, x, y) are renamed to g1, v1, x1, x2
SELECT pp.g1, pp.v1, pp.v2, pp.v3
FROM rasters AS r
CROSS JOIN LATERAL
ROWS FROM (
ST_PixelAsPolygons(rast, 1, exclude_nodata_value => false ),
ST_PixelAsPolygons(rast, 2, exclude_nodata_value => false),
ST_PixelAsPolygons(rast, 3, exclude_nodata_value => false )
) AS pp(g1, v1, x1, y1,
g2, v2, x2, y2,
g3, v3, x3, y3 )
WHERE r.name LIKE '%New York%'
AND ( pp.v1 = 0 OR pp.v2 > 0 OR pp.v3 > 0) ;
Observera
Vi använde CROSS JOIN LATERAL i de här exemplen eftersom vi ville vara tydliga med vad vi gör. Eftersom alla dessa exempel är returnerande funktioner kan du ersätta CROSS JOIN LATERAL med , som en förkortning. Vi kommer att använda ett , i nästa uppsättning exempel
30.6.3. Dumpning av polygoner med ST_DumpAsPolygons¶
Raster introducerar också en ytterligare sammansatt typ som kallas geomval. Tänk på en geomval som en avkomma till en geometri och en raster. Den innehåller en geometri och den innehåller ett pixelvärde.
Du kommer att hitta flera rasterfunktioner som returnerar geomval.
En vanligt förekommande funktion som ger ut geomval är ST_DumpAsPolygons, som returnerar en uppsättning sammanhängande bildpunkter med samma värde som en polygon. Återigen kommer detta som standard endast att kontrollera band 1 och utesluta inga datavärden om du inte åsidosätter. I det här exemplet väljs endast polygoner från band 2. Du kan också tillämpa filter på värdena. För de flesta användningsfall är ST_DumpAsPolygons ett bättre alternativ än ST_PixelAsPolygons eftersom det kommer att returnera mycket färre rader.
Detta ger 6 rader och returnerar polygoner som motsvarar bokstäverna i ”Raster”.
SELECT gv.geom , gv.val
FROM rasters AS r,
ST_DumpAsPolygons(rast, 2) AS gv
WHERE r.name LIKE '%New York%'
AND gv.val = 100;
Observera att den inte returnerar en enda geometri, eftersom den hittar en kontinuerlig uppsättning bildpunkter med samma värde som bildar en polygon. Även om alla dessa värden är desamma är de inte kontinuerliga.
Ett vanligt tillvägagångssätt för att producera mer komplexa geometrier är att gruppera efter värdena och slå ihop dem.
SELECT ST_Union(gv.geom) AS geom , gv.val
FROM rasters AS r,
ST_DumpAsPolygons(rast, 2) AS gv
WHERE r.name LIKE '%New York%'
GROUP BY gv.val;
Detta ger dig 2 rader tillbaka som motsvarar orden ”Raster” och ”Hello”.
30.7. Statistik¶
Det viktigaste att förstå om raster är att de är statistiska verktyg för att lagra data i matriser, som du kanske råkar kunna få att se snygga ut på en skärm.
Du kan hitta en meny med dessa statistiska funktioner i Raster Band Statistics.
30.7.1. ST_SummaryStatsAgg och ST_SummaryStats¶
Om du vill ha all statistik för en uppsättning eller raster kan du använda funktionen ST_SummaryStatsAgg.
Denna fråga tar cirka 10 sekunder och ger dig en sammanfattning av hela tabellen:
SELECT (ST_SummaryStatsAgg(rast, 1, true, 1)).* AS sa
FROM o_3_nyc_dem;
Utmatning:
count | sum | mean | stddev | min | max
-----------+------------+------------------+------------------+-----+-----
246794100 | 4555256024 | 18.4577184948911 | 39.4416860598687 | 0 | 411
(1 row)
Det betyder att vi har många bildpunkter och att vår maxhöjd är 411 fot.
Om du har byggt översikter och bara behöver en grov uppskattning av dina mins, maxs och medelvärden använder du en av dina översikter. Nästa fråga returnerar ungefär samma värden för min, max och medel som den tidigare men på cirka 1 sekund istället för 10.
SELECT (ST_SummaryStatsAgg(rast, 1, true, 1)).* AS sa
FROM o_3_nyc_dem ;
Med den här informationen kan vi nu ställa fler frågor.
30.7.2. ST_Histogram¶
I allmänhet vill du inte ha statistik för hela tabellen, utan istället bara statistik för ett visst område, och i så fall vill du också använda våra gamla vänner ST_Intersects och ST_Clip. Om du också är i behov av en rasterstatistikfunktion som inte har en aggregerad version, vill du ha med dig ST_Union för resan.
För nästa exempel använder vi en annan statistikfunktion ST_Histogram som inte har någon aggregerad motsvarighet, och för denna speciella variant är det en funktion som returnerar en uppsättning. Vi använder samma intresseområde som i några tidigare exempel, men vi måste också använda geometri ST_Transform för att omvandla vår NY State Plane meter geometri till våra NYC State Plane feet rasters. Det är nästan alltid mer performant att transformera geometrin istället för raster och definitivt om din geometri bara är en enda.
SELECT (ST_Quantile( ST_Union( ST_Clip(r.rast, g.geom) ), ARRAY[0.25,0.50,0.75, 1.0] )).*
FROM nyc_dem AS r
INNER JOIN
ST_Transform(
ST_Buffer(ST_Point(586598, 4504816, 26918), 100 ),
2263) AS g(geom)
ON ST_Intersects(r.rast, g.geom);
ovanstående fråga slutförs på mindre än 60 ms och matas ut:
quantile | value
----------+-------
0.25 | 52
0.5 | 57
0.75 | 68
1 | 78
(4 rows)
30.8. Skapa derivativa raster¶
PostGIS raster levereras med ett antal funktioner för att redigera raster. Dessa funktioner används både för att redigera och för att skapa derivata rasterdatauppsättningar. Du hittar dessa listade i Raster Editors och Raster Management.
30.8.1. Transformera raster med ST_Transform¶
De flesta av våra data är i NY State Plane meter (SRID: 26918), men vår DEM-rasterdatauppsättning är i NY State Plane feet (SRID: 2263). För att få ett så smidigt arbetsflöde som möjligt måste våra centrala dataset vara i samma spatiala referenssystem.
Funktionen raster ST_Transform är den funktion som är bäst lämpad för detta jobb.
För att skapa en ny nyc dem-datauppsättning i NY State Plane-mätare gör vi följande:
CREATE TABLE nyc_dem_26918 AS
WITH ref AS (SELECT ST_Transform(rast,26918) AS rast
FROM nyc_dem LIMIT 1)
SELECT r.rid, ST_Transform(r.rast, ref.rast) AS rast, r.filename
FROM nyc_dem AS r, ref;
Ovanstående på mitt system tog cirka 1,5 minuter. För en större datamängd skulle det ta mycket längre tid.
I de tidigare nämnda exemplen användes två varianter av rasterfunktionen ST_Transform. Den första var att få ett referensraster som kommer att användas för att omvandla de andra rasterplattorna för att garantera att alla plattor har samma inriktning. Observera att den andra varianten av ST_Transform som används inte ens tar emot en SRID. Detta beror på att SRID och alla pixelskalor och blockstorlekar läses från referensrastret. Om du använde ST_Transform(rast, srid)-formuläret, skulle alla dina raster kunna komma ut med olika inriktning vilket gör det omöjligt att tillämpa en operation som ST_Union på dem.
Det enda problemet med den ovannämnda ST_Transform-strategin är att när du transformerar, finns det transformerade ofta i andra rutor. Om du tittade tillräckligt noga på ovanstående utdata genom att mata ut rastrernas konvexa skrov, skulle du i nästa exempel se irriterande överlappningar runt gränserna.
SELECT rast::geometry
FROM nyc_dem_26918
ORDER BY rid
LIMIT 100;
som visas i pgAdmin skulle se ut ungefär som:
30.8.2. Använda ST_MakeEmptyCoverage för att skapa jämnt uppdelade raster¶
Ett bättre tillvägagångssätt, om än lite långsammare, är att definiera din egen täckningsstruktur från början med ST_MakeEmptyCoverage och sedan hitta de korsande plattorna för varje ny kakel, och ST_Union dessa och sedan använda ST_Transform(ref, ST_Union…) för att skapa varje kakel.
För detta kommer vi att använda en hel del funktioner, som vi lärt oss om tidigare.
DROP TABLE IF EXISTS nyc_dem_26918;
CREATE TABLE nyc_dem_26918 AS
SELECT ROW_NUMBER() OVER(ORDER BY t.rast::geometry) AS rid,
ST_Union(ST_Clip( ST_Transform( r.rast, t.rast, 'Bilinear' ), t.rast::geometry ), 'MAX') AS rast
FROM (SELECT ST_Transform(
ST_SetSRID(ST_Extent(rast::geometry),2263)
, 26918) AS geom
FROM nyc_dem
) AS g, ST_MakeEmptyCoverage(tilewidth => 256, tileheight => 256,
width => (ST_XMax(g.geom) - ST_XMin(g.geom))::integer,
height => (ST_YMax(g.geom) - ST_YMin(g.geom))::integer,
upperleftx => ST_XMin(g.geom),
upperlefty => ST_YMax(g.geom),
scalex => 3.048,
scaley => -3.048,
skewx => 0., skewy => 0., srid => 26918) AS t(rast)
INNER JOIN nyc_dem AS r
ON ST_Transform(t.rast::geometry, 2263) && r.rast
GROUP BY t.rast;
Upprepa samma övning som tidigare:
SELECT rast::geometry
FROM nyc_dem_26918
ORDER BY rid
LIMIT 100;
sett i pgAdmin har vi inte längre överlappningar:
På mitt system tog detta ~ 10 minuter och returnerade 3879 rader. Efter att tabellen har skapats vill vi göra det vanliga med att lägga till ett spatialt index, primärnyckel och begränsningar enligt följande:
ALTER TABLE nyc_dem_26918
ADD CONSTRAINT pk_nyc_dem_26918 PRIMARY KEY(rid);
CREATE INDEX ix_nyc_dem_26918_st_convexhull_gist
ON nyc_dem_26918 USING gist( ST_ConvexHull(rast) );
SELECT AddRasterConstraints('nyc_dem_26918'::name, 'rast'::name);
ANALYZE nyc_dem_26918;
Vilket bör ta mindre än 2 minuter för detta dataset.
30.8.3. Skapa översiktstabeller med ST_CreateOverview¶
Precis som med vårt ursprungliga dataset skulle det vara användbart att ha översiktstabeller för att påskynda vissa operationer. ST_CreateOverview är en funktion som passar för det ändamålet. Du kan använda ST_CreateOverview för att skapa översikter även om du har glömt att skapa en under raster2pgsql-laddningen eller om du har bestämt dig för att du behöver fler översikter.
Vi skapar översikter på nivå 2 och 3 som vi gjorde med vårt original med hjälp av den här koden.
SELECT ST_CreateOverview('nyc_dem_26918'::regclass, 'rast', 2);
SELECT ST_CreateOverview('nyc_dem_26918'::regclass, 'rast', 3);
Denna process tar tyvärr en stund, och en längre stund ju fler rader du har så ha tålamod. För det här datasetet tog det cirka 3-5 minuter för översiktsfaktor 2 och 1 minut för översiktsfaktor 3.
Funktionen ST_CreateOverView lägger också till de nödvändiga begränsningarna så att kolumnerna visas med full detaljrikedom i katalogerna raster_columns och raster_overviews. Den lägger dock inte till index för dem och lägger inte heller till en rid-kolumn. Rid-kolumnen är förmodligen inte nödvändig om du inte behöver en primärnyckel för att redigera med. Du skulle förmodligen vilja ha ett index som du kan skapa med följande:
CREATE INDEX ix_o_2_nyc_dem_26918_st_convexhull_gist
ON o_2_nyc_dem_26918 USING gist( ST_ConvexHull(rast) );
CREATE INDEX ix_o_3_nyc_dem_26918_st_convexhull_gist
ON o_3_nyc_dem_26918 USING gist( ST_ConvexHull(rast) );
Observera
ST_CreateOverview har ett valfritt argument för att ange urvalsmetoden. Om inget anges används standardmetoden ”NearestNeighbor”, som i allmänhet är snabbast att beräkna men som kanske inte är idealisk. Att ändra samplingsmetoder ligger utanför ramen för denna workshop.
30.9. Skärningspunkten mellan raster och geometrier¶
Det finns ett par funktioner som ofta används för att beräkna korsningspunkter mellan raster och geometrier. Vi har redan sett ST_Clip i aktion, som returnerar korsningspunkten mellan ett raster och en geometri som ett raster, men det finns fler. För punktdata är det vanligaste ST_Value och sedan finns ST_Intersection som har flera överbelastningar, vissa returnerar raster och andra returnerar en uppsättning geomval.
30.9.1. Pixelvärden vid en geometrisk punkt¶
Om du behöver returnera värden från raster baserat på korsningspunkten mellan flera ad hoc-geometripunkter använder du ST_Value eller dess närmaste släkting ST_NearestValue.
SELECT g.geom, ST_Value(r.rast, g.geom) AS elev
FROM nyc_dem_26918 AS r
INNER JOIN
(SELECT id, geom
FROM nyc_homicides
WHERE weapon = 'gun') AS g
ON r.rast && g.geom;
Det här exemplet tar cirka 1 sekund att returnera 2444 rader. Om du använde ST_Intersects i stället för && skulle processen ta cirka 3 sekunder. Anledningen till att ST_Intersects är långsammare är att den utför ytterligare en kontroll, i vissa fall en kontroll pixel för pixel. Om du förväntar dig att alla dina punkter ska representeras med data i din rasteruppsättning och dina raster representerar en täckning (en kontinuerlig uppsättning icke-överlappande rasterplattor), är && i allmänhet ett snabbare alternativ.
Om dina rasterdata inte är tätbefolkade eller om du har överlappande raster (t.ex. de representerar olika observationer i tid), eller om de är skeva (inte axeljusterade), finns det en fördel med att ST_Intersects rensar ut de falska positiva.
30.9.2. ST_Intersection rasterstil¶
Precis som du kan beräkna korsningspunkten mellan två geometrier med hjälp av ST_Intersection, kan du beräkna korsningspunkten mellan två raster eller ett raster och en geometri med hjälp av raster ST_Intersection.
Vad du får ut av det här odjuret är två olika sorters saker:
Skär en geometri med ett raster och du får en uppsättning geomval-avkommor. Kanske en, men oftast många.
Korsa 2 raster och du får ett enda ”raster” tillbaka.
Den gyllene regeln för både rasterkorsning och geometrikorsning är att båda parter måste ha samma spatiala referenssystem. För raster/raster måste de också ha samma inriktning.
Här är ett exempel som svarar på en fråga som du kanske har varit nyfiken på. Om vi skopar våra höjder i 5 hinkar med höjdvärden, vilket höjdintervall resulterar i flest dödsfall i vapen? Vi vet baserat på vår tidigare sammanfattande statistik att 0 är det lägsta värdet och 411 är det högsta värdet för höjd i vår nyc dem dataset, så vi använder det som min och max värde för vår width_bucket samtal.
SELECT ST_Transform(ST_Union(gv.geom),4326) AS geom ,
MIN(gv.val) AS min_elev, MAX(gv.val) AS max_elev,
count(g.id) AS count_guns
FROM nyc_dem_26918 AS r
INNER JOIN nyc_homicides AS g
ON ST_Intersects(r.rast, g.geom)
CROSS JOIN
ST_Intersection( g.geom,
ST_Clip(r.rast,ST_Expand(g.geom, 4) )
) AS gv
WHERE g.weapon = 'gun'
GROUP BY width_bucket(gv.val, 0, 411, 5)
ORDER BY width_bucket(gv.val, 0, 411, 5);
Finns det ett viktigt samband mellan mord med skjutvapen och höjd över havet? Förmodligen inte.
Låt oss ta en titt på raster / raster-korsning:
SELECT ST_Intersection(r1.rast, 1, r2.rast, 1, 'BAND1')
FROM nyc_dem_26918 AS r1
INNER JOIN
rasters AS r2 ON ST_Intersects(r1.rast,1, r2.rast, 1);
Vad vi får är två rader med NULLLs, och om du har din PostgreSQL inställd på att visa meddelanden ser du det:
OBS: De två rastrerna som tillhandahålls har inte samma inriktning. Återlämnar NULL
För att åtgärda detta kan vi anpassa den ena till den andra när den kommer ut ur porten med hjälp av ST_Resample.
SELECT ST_Intersection(r1.rast, 1,
ST_Resample( r2.rast, r1.rast ), 1,
'BAND1')
FROM nyc_dem_26918 AS r1
INNER JOIN
rasters AS r2 ON ST_Intersects(r1.rast,1, r2.rast, 1);
Låt oss också sammanställa det till en enda statistikpost
SELECT (
ST_SummaryStatsAgg(
ST_Intersection(r1.rast, 1,
ST_Resample( r2.rast, r1.rast ), 1, 'BAND1'),
1, true)
).*
FROM nyc_dem_26918 AS r1
INNER JOIN
rasters AS r2 ON ST_Intersects(r1.rast,1, r2.rast, 1);
som matar ut:
count | sum | mean | stddev | min | max
-------+-------+-----------------+------------------+-----+-----
2075 | 99536 | 47.969156626506 | 9.57974836865737 | 33 | 62
(1 row)
30.10. Kartlägga algebrafunktioner¶
Map algebra är tanken att man kan göra matematik på sina pixelvärden. Funktionerna ST_Union och ST_Intersection som nämndes tidigare är ett snabbt specialfall av kartalgebra. Sedan finns det ST_MapAlgebra familjen av funktioner som gör att du kan definiera din egen galna matematik, men på bekostnad av prestanda.
Folk har för vana att hoppa till ST_MapAlgebra, förmodligen för att namnet låter så coolt och sofistikerat. Vem skulle inte vilja berätta det för sina vänner? Jag använder en funktion som heter ST_MapAlgebra Mitt råd är att utforska andra funktioner innan du drar fram hagelgeväret. Ditt liv kommer att bli enklare och din prestanda kommer att bli 100 gånger bättre, och din kod kommer att bli kortare.
Innan vi visar ST_MapAlgebra, utforskar vi andra funktioner som passar in under Map Algebra-familjen av funktioner och har i allmänhet bättre prestanda än ST_MapAlgebra.
30.10.1. Omklassificera ditt raster med hjälp av ST_Reclass¶
En ofta förbisedd map-algebraisk funktion är funktionen ST_Reclass, som sitter i bakgrunden och väntar på att någon ska upptäcka den kraft och snabbhet som den kan erbjuda.
Vad gör ST_Reclass? Som namnet antyder omklassificerar den dina pixelvärden baserat på minimalistisk intervallalgebra.
Låt oss återgå till våra NYC Dems. Kanske bryr vi oss bara om att klassificera våra höjder som 1) låg, 2) medium, 3) hög och 4) riktigt hög. Vi behöver inte 411 värden, vi behöver bara 4. Med det sagt låt oss göra lite omklassificering.
Klassificeringssystemet styrs av reclass-uttrycket.
WITH r AS ( SELECT ST_Union(newrast) As rast
FROM nyc_dem_26918 AS r
INNER JOIN ST_Buffer(ST_Point(586598, 4504816, 26918), 1000 ) AS g(geom)
ON ST_Intersects( r.rast, g.geom )
CROSS JOIN ST_Reclass( ST_Clip(r.rast,g.geom), 1,
'[0-10):1, [10-50):2, [50-100):3,[100-:4','4BUI',0) AS newrast
)
SELECT SUM(ST_Area(gv.geom)::numeric(10,2)) AS area, gv.val
FROM r, ST_DumpAsPolygons(rast) AS gv
GROUP BY gv.val
ORDER BY gv.val;
Vilket skulle ge utdata:
area | val
------------+-----
6754.04 | 1
1753117.51 | 2
1355232.37 | 3
1848.75 | 4
(4 rows)
Om det här var ett klassificeringssystem som vi föredrog skulle vi kunna skapa en ny tabell med ST_Reclass för att räkna om varje bricka.
30.10.2. Färglägga dina raster med ST_ColorMap¶
Funktionen ST_ColorMap är en annan mapalgebraisk funktion som omklassificerar dina pixelvärden. Förutom att det är bandskapande. Det omvandlar ett enda bandraster som våra Dems till ett visuellt presentabelt 3- eller 4-bandat raster.
Du kan använda en av de inbyggda färgkartorna enligt nedan om du inte vill krångla med att skapa en.
SELECT ST_ColorMap( ST_Union(newrast), 'bluered') As rast
FROM nyc_dem_26918 AS r
INNER JOIN
ST_Buffer(
ST_Point(586598, 4504816, 26918), 1000
) AS g(geom)
ON ST_Intersects( r.rast, g.geom)
CROSS JOIN ST_Clip(rast, g.geom) AS newrast;
Vilket ser ut som:
Ju blåare färg desto lägre höjd och ju rödare färg desto högre höjd.
