create table nys_mesonet_raster as with inputs as (
select
-- Sets the pixel size to 0.01 decimal degrees since
-- we are using EPSG 4326
0.01 :: float8 as pixelsize,
-- Sets the smoothing algorithim from "gdal_grid"
'invdist:power:5.5:smoothing:2.0' as algorithm,
st_collect(
-- Creates a geometry collection of our data and forces
-- a Z coodrinate of the max temparature
st_force3dz(
geom,
case
when "apparent_temperature_max [degc]" = '' then null
else "apparent_temperature_max [degc]" :: numeric
end
)
) as geom,
-- Expands the grid to add room aroung the edges
st_expand(st_collect(geom), 0.5) as ext
from
nys_mesonet_202307
where
time_end = '2023-07-23 12:00:00 EDT'
),
sizes AS (
SELECT
ceil((ST_XMax(ext) - ST_XMin(ext)) / pixelsize) :: integer AS width,
ceil((ST_YMax(ext) - ST_YMin(ext)) / pixelsize) :: integer AS height,
ST_XMin(ext) AS upperleftx,
ST_YMax(ext) AS upperlefty
FROM
inputs
)
SELECT
-- Sets 1 as the raster ID since we only have one raster
1 AS rid,
ST_InterpolateRaster(
-- The geometry collection that will be used to interpolate to the raster
geom,
-- The algorithim we defined
algorithm,
ST_SetSRID(
-- This creates the band
ST_AddBand(
-- Creates the empty raster with our arguments from the previous CTEs
ST_MakeEmptyRaster(width, height, upperleftx, upperlefty, pixelsize),
-- Sets the default values as a 32 bit float
'32BF'
),
-- The SRID the raster will use
ST_SRID(geom)
)
) AS rast
FROM
sizes,
inputs