Raster Mashups

[rburghol]

Basic PostGIS raster tools

Combine various raster tools for a date:

• Clip a watersheds best data with st_clip(): https://postgis.net/docs/RT_ST_Clip.html
• Resample all rasters to common resolution with:
  • st_resample(): https://postgis.net/docs/RT_ST_Resample.html
  • or, ST_Rescale(): https://postgis.net/docs/RT_ST_Rescale.html
  • or, st_resize(): https://postgis.net/docs/RT_ST_Resize.html
  • see bug fix using resampling for small watersheds: Bug: issue with NLDAS not generating precip in small/coastal areas #62
• Combine watershed clip rasters with:
  • st_union(): https://postgis.net/docs/RT_ST_Union.html
  • st_mapalgebra(): https://postgis.net/docs/RT_ST_MapAlgebra_expr.html
• and/or Set raster cell values with st_setvalues(): https://postgis.net/docs/RT_ST_SetValues.html
• More info: https://postgis.net/docs/RT_reference.html

Merging Rasters

• ST_UNION() combines rasters, with options for the combining algorithm.
  • The default combining algorithm is to overlay the 1st raster with the values from the 2nd raster.
  • Because st_union is an aggregate function, it is unclear how to specify the order of the input rasters. to be determined.
    • Have also looked at using raster math to combine like so: st_mapalgebra(tpr.rast, 1, ST_Resample(st_clip(w.rast, fgeo.dh_geofield_geom), tpr.rast), 1, '[rast1] + [rast2]') as rast
    • Were rast1 = tpr.rast and rast2 = st_resample(...)

Temporal Disaggregation

• Create 24x1 hour rasters per day with % of daily total P occurring in that hour
• Multiply best guess daily raster by corresponding day hourly rasters to get final model precip

Data Model

• Options:
  • Store as separate variable names for each scenario combination (bas idea)
  • store as model_prop ->

Create Best Fit Raster Dataset

Create Temp Table to import ratings

• note: this could be inserted into the timeseries table as single factors
• This example is for best method file created here: https://github.com/HARPgroup/HARParchive/blob/master/HARP-2024-2025/methods/simple_monthly_lm.R
• Or can be inserted into an amalgamated raster

create temp table week_ratings (
  yr integer,
  mo integer,
  wk integer,
  best_method varchar(16),
  precip_in float8
);

copy week_ratings from '/tmp/met_weekly_ratings_01665500.csv' with HEADER CSV;

• Add the varid column

alter table week_ratings add column best_varid integer;
alter table week_ratings add column best_varkey varchar(64);

update week_ratings set best_varkey = CASE
  WHEN best_method = 'prism' THEN 'prism_mod_daily'
  WHEN best_method = 'daymet' THEN 'daymet_mod_daily'
  ELSE 'do_not_overwrite_nldas2_obs_hourly'
END;

Create a base raster for the model domain

• NLDAS2 resampled to daymet resolution (make smallest possible)
• copy and type and add a band with 0.0 value at each cell

--use an existing raster as template for new raster
-- adds a single band, with value 0.0 in each pixel using the band type as defined by the daymet raster dataset
create table raster_templates(rid serial, varkey varchar(255), rast raster);
INSERT INTO raster_templates(varkey, rast)
SELECT 'daymet_mod_daily', ST_AddBand(ST_MakeEmptyRaster(rast), 1, ST_BandPixelType(rast), 0.0)
FROM dh_timeseries_weather 
WHERE featureid in (
  select hydroid from dh_feature where hydrocode = 'cbp6_met_coverage'
) 
and varid in (select hydroid from dh_variabledefinition where varkey like 'daymet_mod_daily')
and rast is not null 
limit 1;
-- now do nldas2
INSERT INTO raster_templates(varkey, rast)
SELECT 'nldas2_obs_hourly', ST_AddBand(ST_MakeEmptyRaster(rast), 1, ST_BandPixelType(rast), 0.0)
FROM dh_timeseries_weather 
WHERE featureid in (
  select hydroid from dh_feature where hydrocode = 'cbp6_met_coverage'
) 
and varid in (select hydroid from dh_variabledefinition where varkey like 'nldas2_obs_hourly')
and rast is not null 
limit 1;

• Verify Extent of base raster

select st_astext(st_envelope(rast)) from raster_templates where varkey = 'daymet_mod_daily';
-- POLYGON((-83.66607521306385 44.095589132351506,-73.52855441054511 44.095589132351506,-73.52855441054511 33.95806832983277,-83.66607521306385 33.95806832983277,-83.66607521306385 44.095589132351506))

Create a raster clipped to the coverage.

• right now they are just daily
• later we may have a daily/hourly raster for each data source pre-generated and interpolated
• todo: see if we can use st_union() to combine each catchment best-fit raster to overlay the base raster, instead of st_mapalgebra(), since the default behavior of st_union is to simply replace the cells from the first raster with those in the 2nd raster -- this will work, but only if our subsequent rasters are clipped to their watershe boundaries, which I think is currently the case (see the image of the 01665500 raster below).
• this code uses st_resample() to insure matching raster cells, but this could be done in the import workflow, by simply running an SQL update query on the resulting data, and annotating this in the varkey info
  • Omitting this step will throw an error if the two rasters are not exactly same precision ERROR: rt_raster_iterator: The set of rasters provided (custom extent included, if appropriate) do not have the same alignment

create temp table raster_best(tid serial, featureid integer, tstime bigint, tsendtime bigint, varid integer, rast raster);
insert into raster_best(featureid, tstime, tsendtime, varid, rast)
select f.hydroid, w.tstime, w.tsendtime, w.varid, 
  st_mapalgebra(
      tpr.rast, 1, 
      ST_Resample(st_clip(w.rast, fgeo.dh_geofield_geom), tpr.rast), 1, 
      '[rast1] + [rast2]'
    ) as rast
from dh_feature as f
left outer join field_data_dh_geofield as fgeo
on (
  fgeo.entity_id = f.hydroid
  and fgeo.entity_type = 'dh_feature'
)
left outer join week_ratings as r
on (
 yr = 1984
)
left outer join dh_variabledefinition as v
on (
  v.varkey = r.best_varkey
)
left outer join dh_timeseries_weather as w
on (
   extract( year from to_timestamp(w.tsendtime) ) = r.yr
   and extract( month from to_timestamp(w.tsendtime) ) = r.mo
   and extract( week from to_timestamp(w.tsendtime) ) = r.wk 
  and v.hydroid = w.varid
)
left outer join raster_templates as tpr
on ( varkey = 'daymet_mod_daily')
where f.hydrocode = 'usgs_ws_01665500'
and w.tid is not null;

Export Sample Records as geotiff or ping

Using gdal_translate from cmd line

• Must have the rasters in a non-temp table (not a worry, as we will add these to the dh_timeseries_weather table eventually)
  • create temp table raster_best_out as select * from raster_best;
• Now export

PGHOST=dbase2 PGPORT=5432 PGUSER=rob gdal_translate -of GTiff "PG:dbname=drupal.dh03 table=raster_best_out where='tid = 1 ' " outfile.tiff

Using postgis internal native stuff

• This is an interesting piece of code from SO

SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
select tid, lowrite(lo_open(tid, 131072), tiff) As num_bytes
FROM
( VALUES (lo_create(0),
ST_Astiff( (SELECT rast FROM raster_best where tid = 1))
) ) As v(tid,tiff);

• Which returns:

   tid   | num_bytes
---------+-----------
 9155256 |      4002
(1 row)

• which can then be used to save to file (note this file gets saved in the /tmp directory on deq2, even though the database is dbase2):

\lo_export 9155256 '/tmp/myraster.tiff'
SELECT lo_unlink(9155256 );

• copy with cp /tmp/myraster.tiff www/files/met/test_01665500.tiff
• view at: http://deq1.bse.vt.edu:81/files/met/test_01665500.tiff

Image 1: Best fit raster for 01665500, January 1 1984.