Create `seg_maps` using postgis for use with C++ routines

[rburghol]

@gopal-bhatt @COBrogan - I did this not for the tidal segments but because I had located a southern rivers land seg that never got handled properly, in our phase 5 met project. This is a fast and simple way to insure that we have up to date met data raster pixel maps for the C++ style functions (i.e. NLDAS2GRIB ...). It is pretty fast, only a few seconds per land seg to do the intersection, and since it only has to be done once, it is sufficient IMO. It is still experimental, that is, the results that I am getting are somewhat different than the analysts GIS methods from the original project. But, that could be an improvement for all I know, we need to QA (the x values are nearly identical, but the y values are off by 30-40 positions)

• The below is tested for a phase 5 segment that we somehow missed in a project several years ago (or deleted by accident).
  • A51149
  • entry: A51149 12 x382y127 x383y127 x384y127 x385y127 x382y128 x383y128 x384y128 x385y128 x382y129 x383y129 x384y129 x385y129
  • On our modeling system this will work: /opt/model/model_meteorology/sh/segmap_gis A51149 cbp532_landseg

Final Query (see below for table setup)

\set landseg 'A51149'
\set ftype 'cbp532_landseg'
\pset pager 0

select :'landseg' || ' ' || count(seg_map.xy) || ' ' || string_agg(seg_map.xy, ' ')
from (
  select 'x' || x || 'y' || y as xy
  from dh_feature_fielded as f
  left outer join (
    select (ST_PixelAsPoints(rast, 1)).*
    from rastemp
  ) as rastemp
  on (
    dh_geofield_geom && st_setsrid(rastemp.geom, 4326)
  ) 
  where hydrocode = :'landseg' and ftype = :'ftype'
) as seg_map;

Bash script of query

#!/bin/bash
. om_config

if [ $# -lt 2 ]; then
  echo 1>&2 "Use: segmap_gis landseg ftype"
  echo 1>&2 "Ex: segmap_gis A51149 cbp532_landseg"
  exit
fi

sql="
\\set landseg '$1' \n
\\set ftype '$2' \n
\\pset pager 0 \n

select :'landseg' || ' ' || count(seg_map.xy) || ' ' || string_agg(seg_map.xy, ' ')
from (
  select 'x' || x || 'y' || y as xy
  from dh_feature_fielded as f
  left outer join (
    select (ST_PixelAsPoints(rast, 1)).*
    from rastemp
  ) as rastemp
  on (
    dh_geofield_geom && st_setsrid(rastemp.geom, 4326)
  ) 
  where hydrocode = :'landseg' and ftype = :'ftype'
) as seg_map;
set -f
echo -e $sql | psql  $DB_HOST $DB_NAME

Table Setup (really only needs to be done once, but should have a better name than rastemp if it is to be permanent)

# create a full export from nldas2 as sql
raster2pgsql /backup/meteorology/2000/001/NLDAS_FORA0125_H.A20000101.0100.002.grb rastemp > rastemp.sql
# import as table rastemp
cat rastemp.sql | psql -h dbase2 drupal.dh03
# psql
psql -h dbase2 drupal.dh03

• After setting up the raster, find the points

select rastemp.*
from dh_feature_fielded as f
left outer join (
  select (ST_PixelAsPoints(rast, 1)).*
  from rastemp
) as rastemp
on (
  dh_geofield_geom && st_setsrid(rastemp.geom, 4326)
) where hydrocode = 'A51149' and ftype = 'cbp532_landseg';

• Results:

                    geom                    |        val         |  x  |  y
--------------------------------------------+--------------------+-----+-----
 010100000079E92631085853C0F2D24D6210A04240 |  5.650000000000034 | 382 | 127
 010100000079E92631085053C0F2D24D6210A04240 |  5.650000000000034 | 383 | 127
 010100000079E92631084853C0F2D24D6210A04240 |  5.670000000000016 | 384 | 127
 010100000079E92631084053C0F2D24D6210A04240 |  5.650000000000034 | 385 | 127
 010100000079E92631085853C0F2D24D6210904240 | 5.8799999999999955 | 382 | 128
 010100000079E92631085053C0F2D24D6210904240 |  5.730000000000018 | 383 | 128
 010100000079E92631084853C0F2D24D6210904240 | 5.7000000000000455 | 384 | 128
 010100000079E92631084053C0F2D24D6210904240 |   5.78000000000003 | 385 | 128
 010100000079E92631085853C0F2D24D6210804240 |  6.150000000000034 | 382 | 129
 010100000079E92631085053C0F2D24D6210804240 |  5.910000000000025 | 383 | 129
 010100000079E92631084853C0F2D24D6210804240 |  5.920000000000016 | 384 | 129
 010100000079E92631084053C0F2D24D6210804240 |  6.050000000000011 | 385 | 129