Portfolio for Scott Parks: Out of DB Rasters

It's often convenient to work with rasters having all of PostGIS at hand. But, having to load large or many rasters (out-of-db or not) for a one time purpose is not desirable. The following works whether your raster(s) are local or remote, with some caveats.

Don't confuse this with loading out-of-db rasters using raster2pgsql. You would never want to do that with remote rasters, COG or not. This is strictly working with external rasters that are not internal to the database in any manner.

For remote rasters it's highly desirable to have them in Cloud Native GeoTiff (COG) format. Retrieving just a small subset of data is potentially much more efficient than grabbing the whole raster. With a remote and large GeoTiff, you should probably download it first. For local rasters, it's not an issue. For formats like NetCDF, GRIB or HDF, converting them to GeoTiff first is required. You will also need some basic information about that raster, such as, width/height and pixel spatial resolution.

The example below extracts elevation from a .vrt file, which is a virtual spatial format. No data exists in the file, only information about 591 Cloud Optimized GeoTiff (COG) files, which at 10 meter resolution covers a significant portion of the U.S. lower 48 states. Let's look at it with GDAL (version 3.12):

gdal info conusDem10m.vrt

Driver: VRT/Virtual Raster
Files: conusDem10m.vrt
       ./dem10/cogs/n25w081.tif
       ./dem10/cogs/n25w082.tif
       ...
  Size is 637212, 280812
Coordinate System is:
GEOGCRS["NAD83",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101004,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4269]]
Data axis to CRS axis mapping: 2,1
Origin = (-125.000555556292994,50.000555556194911)
Pixel Size = (0.000092592592695,-0.000092592592695)
Corner Coordinates:
Upper Left  (-125.0005556,  50.0005556) (125d 0' 2.00"W, 50d 0' 2.00"N)
Lower Left  (-125.0005556,  23.9994444) (125d 0' 2.00"W, 23d59'58.00"N)
Upper Right ( -65.9994444,  50.0005556) ( 65d59'58.00"W, 50d 0' 2.00"N)
Lower Right ( -65.9994444,  23.9994444) ( 65d59'58.00"W, 23d59'58.00"N)
Center      ( -95.5000000,  37.0000000) ( 95d30' 0.00"W, 37d 0' 0.00"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  Min=-49.798 Max=4099.564 
  NoData Value=-999999
  Metadata:
    STATISTICS_APPROXIMATE=YES
    STATISTICS_MAXIMUM=4099.564453125
    STATISTICS_MEAN=742.26566521538
    STATISTICS_MINIMUM=-49.797824859619
    STATISTICS_STDDEV=736.93088840297
    STATISTICS_VALID_PERCENT=58.96

What we are most concerned with from this output is the projection, which is "EPSG",4269 and the pixel size, which is Pixel Size = (0.000092592592695,-0.000092592592695). It is imperative that all COGS have the same pixel size. With that, we can now write our SQL:

with dem as (
   select
      -- 2000, 2000 is the width and height in pixels. -118.4058, 36.761 is the upper left corner of the temporary raster we want to create
      st_transform(ST_AddBand(
         ST_MakeEmptyRaster(2000, 2000, -118.4058, 36.7612, 0.000092592592695, 0.000092592592695, 0, 0, 4269)
         ,1
         ,'./vrt/conusDem10m.vrt'
         ,NULL::int[]
      ), 4269) as rast
)

select 
   distinct on (elev_m, wkt_geometry)
   st_value(d.rast, p.wkb_geometry)::int as elev_m
   ,st_astext(st_snaptogrid(p.wkb_geometry, .0001)) as wkt_geometry
from dem d
join trail p on st_intersects(p.wkb_geometry, d.rast)

limit 10

The output should look something like:

 elev_m |       wkt_geometry       
--------+--------------------------
   3016 | POINT(-118.3963 36.7712)
   3021 | POINT(-118.3956 36.7715)
   3026 | POINT(-118.3989 36.7713)
   3027 | POINT(-118.3971 36.7712)
   3028 | POINT(-118.3989 36.7714)
   3030 | POINT(-118.3991 36.7715)
   3031 | POINT(-118.3944 36.7717)
   3031 | POINT(-118.3948 36.7717)
   3032 | POINT(-118.3988 36.7714)
   3034 | POINT(-118.3926 36.7711)
(10 rows)

If you're working with data locally, it's practical to use any raster format you like, regardless of file size, .vrt, .tif, whatever.

That's it!

External Database Rasters