-- Section: Construct rasters from geometries
-- out geometry in fixed pixel size --
-- <start id="code_bag_o_rasters" /> --
CREATE SCHEMA ch07;
CREATE TABLE ch07.bag_o_rasters(
   rid serial primary key
 , rast_name text, rast raster);
-- <end id="code_bag_o_rasters" />
-- <start id="code_geom_st_asraster_fixed" /> --
INSERT INTO ch07.bag_o_rasters(rast_name, rast)  
WITH a1 AS (   -- <co id="co_geom_st_asraster_fixed_1" /> --
 SELECT 
  ST_Buffer(ST_GeomFromText('LINESTRING(448252 5414206,448289 5414317
  ,448293 5414330,448324 5414417,448351 5414495)', 32631), 10) 
  As geom ) 
SELECT 'disprop road', ST_AsRaster(geom,50,500)   -- <co id="co_geom_st_asraster_fixed_2" /> --
FROM a1 
UNION ALL
SELECT 'proport fixed w road'
	, ST_AsRaster(geom        -- <co id="co_geom_st_asraster_fixed_3" /> --
	   ,200  -- <co id="co_geom_st_asraster_fixed_4" /> --
	   ,((ST_YMax(geom) - ST_YMin(geom)) *     
	     200/(ST_XMax(geom) - ST_XMin(geom)))::integer  -- <co id="co_geom_st_asraster_fixed_5" /> --
	) 
FROM a1;
-- <end id="code_geom_st_asraster_fixed" /> --
-- [1] road in paris 
-- [2] forced to 50x500 pixel
-- [3] proportional fixed width
-- [4] width 200px 
-- [5] height proportional to width
-- first load with raster2pgsql --
-- <start id="code_stamp_raster_as_reo" /> --
UPDATE ch07.bag_o_rasters SET rast_name = 'reo' 
WHERE rast_name IS NULL;
-- <end id="code_stamp_raster_as_reo" /> --

-- <start id="code_geom_st_asraster_aligned" /> --
INSERT INTO ch07.bag_o_rasters(rast_name, rast)  
WITH ref AS(  -- <co id="co_geom_st_asraster_aligned_1" /> --
	SELECT ST_MakeEmptyRaster( 500, 500, 445000, 5415000, 2, -2
	        , 0, 0, 32631) As rast
	)
, a1 AS (   -- <co id="co_geom_st_asraster_aligned_2" /> --
 SELECT 
  ST_Buffer(ST_GeomFromText('LINESTRING(448252 5414206,448289 5414317
  ,448293 5414330,448324 5414417,448351 5414495)', 32631), 10) 
  As geom ) 
SELECT 'canvas aligned road'
   , ST_AsRaster(geom, ref.rast, '8BUI'::text)   -- <co id="co_geom_st_asraster_aligned_3" /> --
FROM a1 CROSS JOIN ref;
-- <end id="code_geom_st_asraster_aligned" /> --
-- [1] define reference raster 500x500 starting at 445000, 5415000
-- [2] define road
-- [3] create 8BUI raster (represents subsection of reference)



-- <start id="code_rast_scratch" /> --  
INSERT INTO ch07.bag_o_rasters(rast_name, rast)  
SELECT 'raster 1 band scratch'
  , ST_AddBand( -- <co id="co_code_rast_scratch_1" /> --
	 ST_MakeEmptyRaster( 500, 500, 445000, 5415000, 2
	    , -2, 0, 0, 32631)  -- <co id="co_code_rast_scratch_2" /> --
		, '8BUI'::text, 255, 0 )  -- <co id="co_code_rast_scratch_3" /> --
			As rast  ;
-- <end id="code_rast_scratch" /> --
-- [1] add band 
-- [2] empty raster
-- [3] initialized new band 8BUI,100 for pix value

-- <start id="code_rast_st_setvalue_cr" /> --  
UPDATE ch07.bag_o_rasters AS b
	SET rast = ST_SetValue(rast
				,1 -- <co id="co_code_rast_st_setvalue_cr_1" /> --
				,10 -- <co id="co_code_rast_st_setvalue_cr_2" /> --
				,20 -- <co id="co_code_rast_st_setvalue_cr_3" /> --
				,146 -- <co id="co_code_rast_st_setvalue_cr_4" /> --
			)
	WHERE b.rast_name = 'raster 1 band scratch';
-- <end id="code_rast_st_setvalue_cr" /> --
-- [1] band 1
-- [2] column 10
-- [3] row 20
-- [4] pix val: 146

-- <start id="code_rast_heatmap" /> --  
INSERT INTO ch07.bag_o_rasters(rast_name, rast)  
WITH heatmap As ( -- <co id="co_code_rast_heatmap_1" /> --
	SELECT array_agg( -- <co id="co_code_rast_heatmap_2" /> --
		(ST_Buffer(  
			ST_Translate(
				ST_SetSRID(ST_Point(445500, 5414500)
				  ,32631), -500 + i*150
				   ,-200 + 160*i)
					, i*50)  -- <co id="co_code_rast_heatmap_3" /> --
					,  50 + i*15.0 -- <co id="co_code_rast_heatmap_4" /> --
					)::geomval  -- <co id="co_code_rast_heatmap_5" /> --
				) As gvals -- <co id="co_code_rast_heatmap_2c" /> --
	FROM generate_series(-3,4) As i
)
SELECT 'raster 1 band heatmap'
  , ST_SetValues(rast 
  		, 1  -- <co id="co_code_rast_heatmap_6" /> --
  		, heatmap.gvals) As rast -- <co id="co_code_rast_heatmap_7" /> --
FROM ch07.bag_o_rasters As b CROSS JOIN heatmap
WHERE b.rast_name = 'raster 1 band scratch';
-- <end id="code_rast_heatmap" /> --
-- [1] Define a cte heatmap
-- [2] cte has single column gvals collected geomval
-- [3] geoms are buffered translated points
-- [4] vals are 50 + i*15.0
-- [5] cast each geometry/val to type geomval
-- [6] band 1
-- [7] burn the geom values into band

-- crop and get one band --
-- <start id="code_rast_st_band_clip" /> --
INSERT INTO ch07.bag_o_rasters(rast_name, rast)
SELECT 'reo 1 band crop'
	, ST_Clip(      -- <co id="co_code_rast_st_band_clip_1" /> --
		rast , 1  -- <co id="co_code_rast_st_band_clip_2" /> --
        , ST_Buffer(
            ST_Centroid(rast::geometry)
         , 75)  -- <co id="co_code_rast_st_band_clip_3" /> --
         , 255 -- <co id="co_code_rast_st_band_clip_4" />  --
      )
FROM ch07.bag_o_rasters 
WHERE rast_name = 'reo';
-- <end id="code_rast_st_band_clip" /> --
--[1] clip
--[2] band 1
--[3] center circle
--[4] no data value

-- Listing: disoriented reo  --
-- <start id="code_rast_st_band_shuffle" /> --
INSERT INTO ch07.bag_o_rasters(rast_name, rast)
SELECT 'reo band shuffle'
	, ST_Band(  
		rast 
		, ARRAY[3,1,2]) 

FROM ch07.bag_o_rasters 
WHERE rast_name = 'reo';
-- <end id="code_rast_st_band_shuffle" /> --


-- Listing: Raster meta data
-- <start id="code_basic_raster_meta_data" /> --
SELECT rid As r, rast_name
 , ST_Width(rast) As w, ST_Height(rast) As h
 , round(ST_PixelWidth(rast)::numeric,4) AS pw
 , round(ST_PixelHeight(rast)::numeric,4) As ph
 , ST_SRID(rast) AS srid
 , ST_BandPixelType(rast,1) AS bt
FROM ch07.bag_o_rasters;
-- <end id="code_basic_raster_meta_data" />

-- <start id="code_basic_raster_st_metadata_bandmetadata" /> -- 
SELECT rid As r
  , (rm).upperleftx As ux  -- <co id="co_code_basic_raster_st_metadata_bandmetadata_1a" /> --
  , (rm).numbands As nb -- <co id="co_code_basic_raster_st_metadata_bandmetadata_1b" /> --
  , (rbm).*  -- <co id="co_code_basic_raster_st_metadata_bandmetadata_2" /> --
FROM (SELECT rid, ST_MetaData(rast) As rm
	, ST_BandMetaData(rast,1) As rbm
FROM ch07.bag_o_rasters) As r;
-- <end id="code_basic_raster_st_metadata_bandmetadata" /> --
-- [1 a,b] sub properties of raster meta data
-- [2] all properties of raster band meta data

-- <start id="code_st_aspng" /> --
SELECT ST_AsPNG(rast) As png
FROM ch07.bag_o_rasters 
WHERE rast_name = 'raster 1 band heatmap';
-- <end id="code_st_aspng" /> --



-- <start id="code_raster_georeference_processing_effects" /> --
WITH r As 
(SELECT rast  FROM ch07.bag_o_rasters  -- <co id="co_code_raster_georeference_processing_effects_1" /> --
WHERE  rast_name = 'canvas aligned road')
, r2 AS (
 SELECT 'orig' As op
 , ST_MetaData(rast) As rm
 FROM r
 UNION ALL
 SELECT 'resamp' AS op
  , ST_MetaData(ST_Resample(rast,300,300)) As rm -- <co id="co_code_raster_georeference_processing_effects_2" /> --
  FROM r
 UNION ALL
 SELECT 'tform' AS op
  , ST_MetaData(ST_Transform(rast,4326)) As rm -- <co id="co_code_raster_georeference_processing_effects_3" /> --
  FROM r
  UNION ALL
 SELECT 'resize' AS op
  , ST_MetaData(ST_Resize(rast,0.5,0.5)) As rm -- <co id="co_code_raster_georeference_processing_effects_4" /> --
  FROM r
  UNION ALL
 SELECT 'rescale' AS op
  , ST_MetaData(ST_Rescale(rast,0.5,-0.5)) As rm -- <co id="co_code_raster_georeference_processing_effects_5" /> --
  FROM r
  ) 
SELECT op 
  ,(rm).srid 
  , (rm).width::text || 'x' || (rm).height::Text as wh
  ,(rm).scalex::numeric(7,5)::text || ',' 
    || (rm).scaley::numeric(7,5)::text as sxy
  , (rm).upperleftx::numeric(11,2)::text 
    || ',' || (rm).upperlefty::numeric(12,2)::text As uplxy -- <co id="co_code_raster_georeference_processing_effects_6" /> --
FROM r2;
-- <end id="code_raster_georeference_processing_effects" /> --
--[1] reference raster
--[2] resample 300x300 pixels
--[3] transform to long lat
--[4] resize to 50%
--[5] rescale to 0.5, -0.5 meters/pix
--[6] upper left corner 

-- Listing writeable formats --
-- <start id="code_gdaldrivers_writeable_formats" /> --
SELECT  short_name, long_name 
from ST_GDALDrivers()
ORDER BY short_name;
-- <end id="code_gdaldrivers_writeable_formats" /> --

-- options for particular raster type --
-- <start id="code_gdaldrivers_writeable_formats_options" /> --
SELECT (xpath('@name', g.opt))[1]::text As oname,
       (xpath('@type', g.opt))[1]::text As otype,
       (xpath('@description', g.opt))[1]::text As descrip
FROM (SELECT unnest(
	xpath('/CreationOptionList/Option'
	 , create_options::xml)) As opt
FROM  ST_GDALDrivers()
WHERE short_name = 'USGSDEM') As g;
-- <end id="code_gdaldrivers_writeable_formats_options" /> --

-- write using ST_AsGDALRaster --
-- <start id="code_st_asgdalraster_usgsdem" /> --
SELECT ST_AsGDALRaster(
	ST_Band(rast,1) -- <co id="co_code_st_asgdalraster_usgsdem_1" />  --
  , 'USGSDEM'  -- <co id="co_code_st_asgdalraster_usgsdem_2" />  --
  , ARRAY[ -- <co id="co_code_st_asgdalraster_usgsdem_3" />  --
  'PRODUCER=' || quote_literal('postgis_in_action') 
   , 'INTERNALNAME=' || quote_literal(rast_name)]
   ) As dem
FROM ch07.bag_o_rasters 
WHERE rast_name='raster 1 band heatmap';
-- <end id="code_st_asgdalraster_usgsdem" /> --
--[1] band 1
--[2] format: USDEM
--[3] metadata to stamp in file

-- output with psql --
-- <start id="code_st_asgdalraster_output_psql" /> --
SELECT oid, lowrite(lo_open(oid, 131072), img) As num_bytes
 FROM 
 ( VALUES (lo_create(0) -- <co id="co_code_st_asgdalraster_output_psql_1a" /> --
  ,  (   -- <co id="co_code_st_asgdalraster_output_psql_2a" /> --
   SELECT ST_AsGDALRaster(
	ST_Band(rast,1) 
  , 'USGSDEM' 
  , ARRAY['PRODUCER=' || quote_literal('postgis_in_action') 
   , 'INTERNALNAME=' || quote_literal(rast_name)]
   ) As dem
FROM ch07.bag_o_rasters 
WHERE rast_name='raster 1 band heatmap'
)   -- <co id="co_code_st_asgdalraster_output_psql_2b" /> --
  ) ) As v(oid,img);
  
\lo_export 9585208 'C:/temp/heatmap.dem'    -- <co id="co_code_st_asgdalraster_output_psql_3" /> --
SELECT lo_unlink(9585208);   -- <co id="co_code_st_asgdalraster_output_psql_4" /> --
;
-- <end id="code_st_asgdalraster_output_psql" /> --
--[1] create large object
--[2] query output to dem
--[3] use oid returned by previous query
--[4] when done delete large object


-- Listing histogram breakout --
-- <start id="code_st_histogram_single" /> --
SELECT  (stats).*
FROM (SELECT ST_Histogram(rast
	, 2 -- <co id="co_code_st_histogram_single_1" />  --
	,6 -- <co id="co_code_st_histogram_single_2" />  --
	) As stats
    FROM ch07.bag_o_rasters WHERE rast_name = 'reo') As foo;
-- <end id="code_st_histogram_single" /> --
-- [1] band 2
-- [2] 6 bins

-- informal summary stats --
-- <start id="code_st_summarystats_single" /> --
SELECT  (stats).*
FROM (SELECT ST_SummaryStats(rast
	, 2 
	) As stats
    FROM ch07.bag_o_rasters WHERE rast_name = 'reo') As foo;
-- <end id="code_st_summarystats_single" /> --

-- informal georeferencing --
-- <start id="code_st_georeference_single" /> --
UPDATE ch07.bag_o_rasters 
  SET rast = ST_SetSRID(
  	ST_SetGeoReference(rast, '1 0 0 -1 445139 5415000')
  	,32631)
WHERE rast_name = 'reo 1 band crop';
-- <end id="code_st_georeference_single" /> --


-- Listing code_st_reclass_band2_reo --
-- <start id="code_st_reclass_band2_reo" /> --
INSERT INTO ch07.bag_o_rasters(rast_name, rast)
SELECT 'reo 1 banded band 2 reclass', 
	ST_Reclass(
		ST_Band(rast,2) -- <co id="co_code_st_reclass_band2_reo_1" /> --
	,1 -- <co id="co_code_st_reclass_band2_reo_2" /> --
	, '[0-66]:0,(66-255]:255'::text -- <co id="co_code_st_reclass_band2_reo_3" /> --
	, '8BUI'::text
	, 255 -- <co id="co_code_st_reclass_band2_reo_4" /> --
  ) 
from ch07.bag_o_rasters
WHERE rast_name = 'reo';
-- <end id="code_st_reclass_band2_reo" /> --
-- [1] create new rast with just band 2
-- [2] band 1 of single band
-- [3] reclass so val > 66 and val <= 255 are set to no data
-- [4] set 255 to no data

-- Compare original reo with reclassified --
-- <start id="code_st_convexhull_reo" /> --
SELECT rast_name::varchar(10)
    , ST_AsText(ST_ConvexHull(rast)) As hull
    , ST_AsText(ST_MinConvexHull(rast)) As minhull
FROM ch07.bag_o_rasters
WHERE rast_name IN('reo'
 , 'reo 1 banded band 2 reclass');
-- <end id="code_st_convexhull_reo" /> --
