Changes between Version 77 and Version 78 of PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools

Timestamp:

Jul 28, 2012, 12:24:06 AM (12 years ago)

Author:

qliu

Comment:

—

PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools

@@ -56 +56 @@
 <LI><A href=http://trac.osgeo.org/postgis/wiki/PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools#Week7Report>Week 7 Report (July 2nd - July 6th)</A></LI>
 <LI><A href=http://trac.osgeo.org/postgis/wiki/PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools#Week8Report>Week 8 Report (July 9th - July 13th)</A></LI>
+<LI><A href=http://trac.osgeo.org/postgis/wiki/PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools#Week10Report>Week 10 Report (July 23th - July 27th)</A></LI>
 </UL>
 <br>
@@ -833 +834 @@
 
  * I will be taking a vacation next week. And I will have to speed it up and continue to work on implementing a plpgsql prototype for Approach 3 in the document after I come back.
+
+
+----
+{{{
+#!html
+<div style='float: right; width:100px;'>
+<a href="#PostGISRasterSoCIdea2012-DistanceAnalysisToolsforPostGISRaster">&nbsp;&nbsp;&nbsp;<img src="https://lh3.googleusercontent.com/-YR7gUJULQrM/T8nEPQ-AtzI/AAAAAAAAA1M/T_6vRjE8E4Q/s40/scroll_to_top_icon_40x40.png"><br>back to top</a>
+</div>
+}}}
+
+== Week 10 Report ==
+
+'''What did you get done this week?'''
+ * Wrote plpgsql prototype for implementing Euclidean distance with KNN indexing approach.
+{{{
+#!sql
+----------------------------------------------------------------------
+-- $Id: st_euclideandistance.sql 2012-07-28 $
+----------------------------------------------------------------------
+
+--------------------------------------------------------------------
+-- ST_EuclideanDistance - (for one raster) Return a raster generated from a one band reference raster
+--                        which values are the Euclidean distance
+--                        from pixel to the points (for now only concerning point geometries) 
+--                        from the input source table of geometries 
+--                        (points for now, lines or polygons in the future).
+-- Arguments 
+-- refrast raster -  Reference raster align with which Euclidean distance raster will be created.
+-- pixeltype text - Pixeltype assigned to the resulting raster. By default 
+--                  to be the pixeltype of the reference raster if specified.
+-- sourceschema text - Database schema of the source geometry table.
+-- sourcetable text - Source geometry table name.
+-- sourcegeomcolumn text - Source geometry table geometry column name.
+-- maxdist double precision - Maximum distance from pixel to the source
+--                                when the source is too far from the pixel.
+--                                Pixel that exceeds it will be assigned nodatavalue.
+--------------------------------------------------------------------
+DROP FUNCTION IF EXISTS ST_EuclideanDistance(refrast raster, pixeltype text, sourceschema text, sourcetable text, sourcegeomcolumn text, maxdist double precision);
+CREATE OR REPLACE FUNCTION ST_EuclideanDistance(refrast raster, pixeltype text, sourceschema text, sourcetable text, sourcegeomcolumn text, maxdist double precision) 
+    RETURNS raster AS 
+    $$
+    DECLARE
+        width integer;
+        height integer;
+        x integer;
+        y integer;
+                sourcex float8;
+                sourcey float8;
+                sourcexr integer;
+                sourceyr integer;
+                newx float8;
+                newy float8;
+                newsrid integer;
+                exesql text;
+                sourcegeom geometry;
+        newnodatavalue float8;
+                newinitialvalue float8;
+        newrast raster;
+        newval float8;
+        newpixeltype text;
+        newhasnodatavalue boolean := FALSE;
+                -- Assuming reference raster has only one band
+                band integer := 1;
+    BEGIN
+        -- Check if reference raster is NULL
+        IF refrast IS NULL THEN
+            RAISE NOTICE 'ST_EuclideanDistance: Reference raster is NULL. Returning NULL';
+            RETURN NULL;
+        END IF;
+
+        width := ST_Width(refrast);
+        height := ST_Height(refrast);
+                newsrid := ST_SRID(refrast);
+
+        -- Create a new empty raster having the same georeference as the provided raster
+        newrast := ST_MakeEmptyRaster(width, height, ST_UpperLeftX(refrast), ST_UpperLeftY(refrast), ST_ScaleX(refrast), ST_ScaleY(refrast), ST_SkewX(refrast), ST_SkewY(refrast), newsrid);
+
+        -- Check if the new raster is empty (width = 0 OR height = 0)
+        IF ST_IsEmpty(newrast) THEN 
+            RAISE NOTICE 'ST_EuclideanDistance: Raster is empty. Returning an empty raster';
+            RETURN newrast;
+        END IF;
+        
+        -- Set the new pixeltype
+        newpixeltype := pixeltype;
+                -- If pixeltype not specified then use the pixeltype of the reference raster
+        IF newpixeltype IS NULL THEN
+            newpixeltype := ST_BandPixelType(refrast, band);
+        ELSIF newpixeltype != '1BB' AND newpixeltype != '2BUI' AND newpixeltype != '4BUI' AND newpixeltype != '8BSI' AND newpixeltype != '8BUI' AND 
+               newpixeltype != '16BSI' AND newpixeltype != '16BUI' AND newpixeltype != '32BSI' AND newpixeltype != '32BUI' AND newpixeltype != '32BF' AND newpixeltype != '64BF' THEN
+            RAISE EXCEPTION 'ST_EuclideanDistance: Invalid pixeltype "%". Aborting.', newpixeltype;
+        END IF;
+
+        -- Check for notada value
+        newnodatavalue := ST_BandNodataValue(refrast, band);
+        IF newnodatavalue IS NULL OR newnodatavalue < ST_MinPossibleValue(newpixeltype) OR newnodatavalue > (-ST_MinPossibleValue(newpixeltype) - 1) THEN
+            RAISE NOTICE 'ST_EuclideanDistance: Reference raster do not have a nodata value or is out of range for the new raster pixeltype, nodata value for new raster set to the min value possible';
+            newnodatavalue := ST_MinPossibleValue(newpixeltype);
+        END IF;
+        -- We set the initial value of the resulting raster to nodata value. 
+        -- If nodatavalue is null then the raster will be initialise to ST_MinPossibleValue 
+        newinitialvalue := newnodatavalue;
+        
+        -- Create the raster receiving all the computed distance values. Initialize it to the new initial value.
+        newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue);
+
+                -- Scan pixels in the raster set pixel value to zero if pixel is source point
+                -- There is place for optimization here for a more efficient scanning method
+                FOR x IN 1..width LOOP
+                        FOR y IN 1..height LOOP
+                                exesql := "SELECT " || sourcegeomcolumn || " FROM " || sourcetable || ";";
+                                FOR geom IN EXECUTE(exesql) LOOP
+                                        sourcex := ST_X(geom);
+                                        sourcey := ST_Y(geom);
+                                        sourcexr := ST_World2RasterCoordX(newrast, sourcex, sourcey);
+                                        sourceyr := ST_World2RasterCoordY(newrast, sourcex, sourcey);
+                                        IF x = sourcexr AND y = sourceyr THEN
+                                                newrast := ST_SetValue(newrast, band, x, y, 0);
+                                        END IF;
+                                END LOOP;
+                        END LOOP;
+                END LOOP;
+                
+                -- There is place for optimization here for a more efficient scanning method
+        FOR x IN 1..width LOOP
+            FOR y IN 1..height LOOP
+                                newx := ST_Raster2WorldCoordX(newrast, x, y);
+                                newy := ST_Raster2WorldCoordY(newrast, x, y); 
+                                -- If pixel is not the source point then compute Eculidean distance
+                                IF ST_Value(newrast, band, x, y) <> 0 THEN
+                                        exesql := "SELECT ST_Distance(ST_GeomFromText('POINT(" || newx || " " || newy || ")'," || newsrid || "),ST_Transform(" || sourcegeomcolumn || "," || newsrid || ")) FROM " || sourcetable || " ORDER BY $1 <-> $2 LIMIT 1;";
+                                        newval := EXECUTE(exesql);
+                                END IF;
+                                newrast := ST_SetValue(newrast, band, x, y, newval);
+            END LOOP;
+        END LOOP;
+                
+                -- There is place for optimization here for a more efficient scanning method
+        FOR x IN 1..width LOOP
+            FOR y IN 1..height LOOP
+                                newx := ST_Raster2WorldCoordX(newrast, x, y);
+                                newy := ST_Raster2WorldCoordY(newrast, x, y); 
+                                -- If pixel is the source point then Eculidean distance is zero
+                                IF THEN
+                                        newval := 0;
+                                ELSE
+                                        exesql := "SELECT ST_Distance(ST_GeomFromText('POINT(" || newx || " " || newy || ")'," || newsrid || "),ST_Transform(" || sourcegeomcolumn || "," || newsrid || ")) FROM " || sourcetable || " ORDER BY $1 <-> $2 LIMIT 1;";
+                                        newval := EXECUTE(exesql);
+                                END IF;
+                                newrast := ST_SetValue(newrast, band, x, y, newval);
+            END LOOP;
+        END LOOP;
+        RETURN newrast;
+    END;
+    $$
+    LANGUAGE 'plpgsql';
+}}}
+
+[[BR]][[BR]]
+''' What do you plan on doing next week?'''
+
+ * Test the prototype.
+ * Start writing document for cost-weighted distance.
+ * Start writing C implementation for Euclidean distance.