| | 836 | |
| | 837 | |
| | 838 | ---- |
| | 839 | {{{ |
| | 840 | #!html |
| | 841 | <div style='float: right; width:100px;'> |
| | 842 | <a href="#PostGISRasterSoCIdea2012-DistanceAnalysisToolsforPostGISRaster"> <img src="https://lh3.googleusercontent.com/-YR7gUJULQrM/T8nEPQ-AtzI/AAAAAAAAA1M/T_6vRjE8E4Q/s40/scroll_to_top_icon_40x40.png"><br>back to top</a> |
| | 843 | </div> |
| | 844 | }}} |
| | 845 | |
| | 846 | == Week 10 Report == |
| | 847 | |
| | 848 | '''What did you get done this week?''' |
| | 849 | * Wrote plpgsql prototype for implementing Euclidean distance with KNN indexing approach. |
| | 850 | {{{ |
| | 851 | #!sql |
| | 852 | ---------------------------------------------------------------------- |
| | 853 | -- $Id: st_euclideandistance.sql 2012-07-28 $ |
| | 854 | ---------------------------------------------------------------------- |
| | 855 | |
| | 856 | -------------------------------------------------------------------- |
| | 857 | -- ST_EuclideanDistance - (for one raster) Return a raster generated from a one band reference raster |
| | 858 | -- which values are the Euclidean distance |
| | 859 | -- from pixel to the points (for now only concerning point geometries) |
| | 860 | -- from the input source table of geometries |
| | 861 | -- (points for now, lines or polygons in the future). |
| | 862 | -- Arguments |
| | 863 | -- refrast raster - Reference raster align with which Euclidean distance raster will be created. |
| | 864 | -- pixeltype text - Pixeltype assigned to the resulting raster. By default |
| | 865 | -- to be the pixeltype of the reference raster if specified. |
| | 866 | -- sourceschema text - Database schema of the source geometry table. |
| | 867 | -- sourcetable text - Source geometry table name. |
| | 868 | -- sourcegeomcolumn text - Source geometry table geometry column name. |
| | 869 | -- maxdist double precision - Maximum distance from pixel to the source |
| | 870 | -- when the source is too far from the pixel. |
| | 871 | -- Pixel that exceeds it will be assigned nodatavalue. |
| | 872 | -------------------------------------------------------------------- |
| | 873 | DROP FUNCTION IF EXISTS ST_EuclideanDistance(refrast raster, pixeltype text, sourceschema text, sourcetable text, sourcegeomcolumn text, maxdist double precision); |
| | 874 | CREATE OR REPLACE FUNCTION ST_EuclideanDistance(refrast raster, pixeltype text, sourceschema text, sourcetable text, sourcegeomcolumn text, maxdist double precision) |
| | 875 | RETURNS raster AS |
| | 876 | $$ |
| | 877 | DECLARE |
| | 878 | width integer; |
| | 879 | height integer; |
| | 880 | x integer; |
| | 881 | y integer; |
| | 882 | sourcex float8; |
| | 883 | sourcey float8; |
| | 884 | sourcexr integer; |
| | 885 | sourceyr integer; |
| | 886 | newx float8; |
| | 887 | newy float8; |
| | 888 | newsrid integer; |
| | 889 | exesql text; |
| | 890 | sourcegeom geometry; |
| | 891 | newnodatavalue float8; |
| | 892 | newinitialvalue float8; |
| | 893 | newrast raster; |
| | 894 | newval float8; |
| | 895 | newpixeltype text; |
| | 896 | newhasnodatavalue boolean := FALSE; |
| | 897 | -- Assuming reference raster has only one band |
| | 898 | band integer := 1; |
| | 899 | BEGIN |
| | 900 | -- Check if reference raster is NULL |
| | 901 | IF refrast IS NULL THEN |
| | 902 | RAISE NOTICE 'ST_EuclideanDistance: Reference raster is NULL. Returning NULL'; |
| | 903 | RETURN NULL; |
| | 904 | END IF; |
| | 905 | |
| | 906 | width := ST_Width(refrast); |
| | 907 | height := ST_Height(refrast); |
| | 908 | newsrid := ST_SRID(refrast); |
| | 909 | |
| | 910 | -- Create a new empty raster having the same georeference as the provided raster |
| | 911 | newrast := ST_MakeEmptyRaster(width, height, ST_UpperLeftX(refrast), ST_UpperLeftY(refrast), ST_ScaleX(refrast), ST_ScaleY(refrast), ST_SkewX(refrast), ST_SkewY(refrast), newsrid); |
| | 912 | |
| | 913 | -- Check if the new raster is empty (width = 0 OR height = 0) |
| | 914 | IF ST_IsEmpty(newrast) THEN |
| | 915 | RAISE NOTICE 'ST_EuclideanDistance: Raster is empty. Returning an empty raster'; |
| | 916 | RETURN newrast; |
| | 917 | END IF; |
| | 918 | |
| | 919 | -- Set the new pixeltype |
| | 920 | newpixeltype := pixeltype; |
| | 921 | -- If pixeltype not specified then use the pixeltype of the reference raster |
| | 922 | IF newpixeltype IS NULL THEN |
| | 923 | newpixeltype := ST_BandPixelType(refrast, band); |
| | 924 | ELSIF newpixeltype != '1BB' AND newpixeltype != '2BUI' AND newpixeltype != '4BUI' AND newpixeltype != '8BSI' AND newpixeltype != '8BUI' AND |
| | 925 | newpixeltype != '16BSI' AND newpixeltype != '16BUI' AND newpixeltype != '32BSI' AND newpixeltype != '32BUI' AND newpixeltype != '32BF' AND newpixeltype != '64BF' THEN |
| | 926 | RAISE EXCEPTION 'ST_EuclideanDistance: Invalid pixeltype "%". Aborting.', newpixeltype; |
| | 927 | END IF; |
| | 928 | |
| | 929 | -- Check for notada value |
| | 930 | newnodatavalue := ST_BandNodataValue(refrast, band); |
| | 931 | IF newnodatavalue IS NULL OR newnodatavalue < ST_MinPossibleValue(newpixeltype) OR newnodatavalue > (-ST_MinPossibleValue(newpixeltype) - 1) THEN |
| | 932 | 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'; |
| | 933 | newnodatavalue := ST_MinPossibleValue(newpixeltype); |
| | 934 | END IF; |
| | 935 | -- We set the initial value of the resulting raster to nodata value. |
| | 936 | -- If nodatavalue is null then the raster will be initialise to ST_MinPossibleValue |
| | 937 | newinitialvalue := newnodatavalue; |
| | 938 | |
| | 939 | -- Create the raster receiving all the computed distance values. Initialize it to the new initial value. |
| | 940 | newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue); |
| | 941 | |
| | 942 | -- Scan pixels in the raster set pixel value to zero if pixel is source point |
| | 943 | -- There is place for optimization here for a more efficient scanning method |
| | 944 | FOR x IN 1..width LOOP |
| | 945 | FOR y IN 1..height LOOP |
| | 946 | exesql := "SELECT " || sourcegeomcolumn || " FROM " || sourcetable || ";"; |
| | 947 | FOR geom IN EXECUTE(exesql) LOOP |
| | 948 | sourcex := ST_X(geom); |
| | 949 | sourcey := ST_Y(geom); |
| | 950 | sourcexr := ST_World2RasterCoordX(newrast, sourcex, sourcey); |
| | 951 | sourceyr := ST_World2RasterCoordY(newrast, sourcex, sourcey); |
| | 952 | IF x = sourcexr AND y = sourceyr THEN |
| | 953 | newrast := ST_SetValue(newrast, band, x, y, 0); |
| | 954 | END IF; |
| | 955 | END LOOP; |
| | 956 | END LOOP; |
| | 957 | END LOOP; |
| | 958 | |
| | 959 | -- There is place for optimization here for a more efficient scanning method |
| | 960 | FOR x IN 1..width LOOP |
| | 961 | FOR y IN 1..height LOOP |
| | 962 | newx := ST_Raster2WorldCoordX(newrast, x, y); |
| | 963 | newy := ST_Raster2WorldCoordY(newrast, x, y); |
| | 964 | -- If pixel is not the source point then compute Eculidean distance |
| | 965 | IF ST_Value(newrast, band, x, y) <> 0 THEN |
| | 966 | exesql := "SELECT ST_Distance(ST_GeomFromText('POINT(" || newx || " " || newy || ")'," || newsrid || "),ST_Transform(" || sourcegeomcolumn || "," || newsrid || ")) FROM " || sourcetable || " ORDER BY $1 <-> $2 LIMIT 1;"; |
| | 967 | newval := EXECUTE(exesql); |
| | 968 | END IF; |
| | 969 | newrast := ST_SetValue(newrast, band, x, y, newval); |
| | 970 | END LOOP; |
| | 971 | END LOOP; |
| | 972 | |
| | 973 | -- There is place for optimization here for a more efficient scanning method |
| | 974 | FOR x IN 1..width LOOP |
| | 975 | FOR y IN 1..height LOOP |
| | 976 | newx := ST_Raster2WorldCoordX(newrast, x, y); |
| | 977 | newy := ST_Raster2WorldCoordY(newrast, x, y); |
| | 978 | -- If pixel is the source point then Eculidean distance is zero |
| | 979 | IF THEN |
| | 980 | newval := 0; |
| | 981 | ELSE |
| | 982 | exesql := "SELECT ST_Distance(ST_GeomFromText('POINT(" || newx || " " || newy || ")'," || newsrid || "),ST_Transform(" || sourcegeomcolumn || "," || newsrid || ")) FROM " || sourcetable || " ORDER BY $1 <-> $2 LIMIT 1;"; |
| | 983 | newval := EXECUTE(exesql); |
| | 984 | END IF; |
| | 985 | newrast := ST_SetValue(newrast, band, x, y, newval); |
| | 986 | END LOOP; |
| | 987 | END LOOP; |
| | 988 | RETURN newrast; |
| | 989 | END; |
| | 990 | $$ |
| | 991 | LANGUAGE 'plpgsql'; |
| | 992 | }}} |
| | 993 | |
| | 994 | [[BR]][[BR]] |
| | 995 | ''' What do you plan on doing next week?''' |
| | 996 | |
| | 997 | * Test the prototype. |
| | 998 | * Start writing document for cost-weighted distance. |
| | 999 | * Start writing C implementation for Euclidean distance. |
