PostGIS

Set up PostGIS in database

CREATE EXTENSION IF NOT EXISTS fuzzystrmatch;
CREATE EXTENSION IF NOT EXISTS postgis;
CREATE EXTENSION IF NOT EXISTS postgis_topology;
CREATE EXTENSION IF NOT EXISTS postgis_tiger_geocoder;

SELECT name, default_version,installed_version FROM pg_available_extensions WHERE name LIKE 'postgis%';

PostGIS function raster interpolate value

CREATE OR REPLACE FUNCTION st_ext_valueinterpolated(rast raster, band integer, pt geometry, exclude_nodata_value boolean DEFAULT TRUE) RETURNS float8 AS $$
    DECLARE
        x float8;
        y float8;
        gtype text;
        h_cellsize_x float8;
        h_cellsize_y float8;
        fx float8;
        fy float8;
        ix1 float8;
        iy1 float8;
        ix2 float8;
        iy2 float8;
        nx float8;
        ny float8;
        result float8;

        scaleX float8;
        scaleY float8;
        upperLeftX float8;
        upperLefty float8;

        dx1 float8;
        dx2 float8;
        dy1 float8;
        dy2 float8;
        div float8;
    BEGIN
        gtype := ST_geometrytype(pt);
        IF ( gtype != 'ST_Point' ) THEN
            RAISE EXCEPTION 'Attempting to get the value of a pixel with a non-point geometry';
        END IF;

	IF ST_SRID(pt) != ST_SRID(rast) THEN
            RAISE EXCEPTION 'Raster and geometry do not have the same SRID';
	END IF;
	
	nx := ST_Width(rast);
	ny := ST_Height(rast);
	-- RAISE NOTICE 'i want to print % and %', nx,ny;
        x := ST_x(pt);
        y := ST_y(pt);

	scaleX := ST_ScaleX(rast);
	scaleY := ST_ScaleY(rast);

	upperLeftX := ST_UpperLeftX(rast);
        upperLeftY := ST_UpperLeftY(rast);
	
	h_cellsize_x := scaleX/2;
        h_cellsize_y := scaleY/2;
	
	fx := (x - (upperLeftX + h_cellsize_x))/scaleX;
	fy := (y - (upperLeftY + h_cellsize_y))/scaleY;

	ix1 := floor(fx);
	iy1 := floor(fy);
	
	IF abs(fx - (nx - 1)) < 0.00001 THEN
	  ix1 := ix1 - 1;
	END IF;
	
	IF abs(fy - (ny - 1))  < 0.00001 THEN
	  iy1 :=  iy1- 1;
	END IF;
	  
	ix2 := ix1 + 1;
	iy2 := iy1 + 1;

	IF (ix1 < 0) OR (iy1 < 0) OR (ix2 > nx - 1) OR (iy2 > ny - 1) THEN
		result := ST_BandNoDataValue(rast,band);
    ELSE
		dx1 = x - (upperLeftX + ix1*scaleX + h_cellsize_x);
		dy1 = y - (upperLeftY + iy1*scaleY + h_cellsize_y);
		dx2 = (upperLeftX + ix2*scaleX + h_cellsize_x) - x;
		dy2 = (upperLeftY + iy2*scaleY + h_cellsize_y) - y;
		div = scaleX*scaleY;
	  
		--RAISE NOTICE 'i want to print % ; % --- % ; % --- div %', dx1, dy1, dx2, dy2, div;
	  
		-- offset all values by + 1 as the index in postgres starts at 1,1 instead of 0,0
		ix1 := ix1 + 1;
        iy1 := iy1 + 1;
        ix2 := ix2 + 1;
        iy2 := iy2 + 1;
	  
		--RAISE NOTICE 'i want to print % ; % --- % ; %', ST_Value(rast, ix1::int, iy1::int), ST_Value(rast, ix2::int, iy1::int), ST_Value(rast, ix1::int, iy2::int), ST_Value(rast, ix2::int, iy2::int);
		result := (ST_Value(rast, ix1::int, iy1::int)*dx2*dy2)/div + (ST_Value(rast, ix2::int, iy1::int)*dx1*dy2)/div + (ST_Value(rast, ix1::int, iy2::int)*dx2*dy1)/div + (ST_Value(rast, ix2::int, iy2::int)*dx1*dy1)/div;
    END IF;
	
        RETURN result;
    END;
 $$
 LANGUAGE 'plpgsql' IMMUTABLE STRICT PARALLEL SAFE;