| | 1043 | |
| | 1044 | '''PL/pgSQL script for Euclidean Distance Raster:''' |
| | 1045 | |
| | 1046 | {{{ |
| | 1047 | #!sql |
| | 1048 | |
| | 1049 | ---------------------------------------------------------------------- |
| | 1050 | -- $Id: st_euclideandistance.sql 2012-07-28 $ |
| | 1051 | ---------------------------------------------------------------------- |
| | 1052 | -------------------------------------------------------------------- |
| | 1053 | -- ST_EuclideanDistance - (for one raster) Return a raster generated from a one band reference raster |
| | 1054 | -- which values are the Euclidean distance |
| | 1055 | -- from pixel to the points (for now only concerning point geometries) |
| | 1056 | -- from the input source table of geometries |
| | 1057 | -- (points for now, lines or polygons in the future). |
| | 1058 | -- Arguments |
| | 1059 | -- width,height,rastulx,rastuly,rastscalex,rastscaley,rastskewx,rastskewy |
| | 1060 | -- - Dimenssions and georeference of resulted raster. |
| | 1061 | -- rastsrid - SRID of the resulted raster |
| | 1062 | -- rastndv - Nodata value for resulted raster |
| | 1063 | -- pixeltype text - Pixeltype assigned to the resulting raster. By default |
| | 1064 | -- to be the pixeltype of the reference raster if specified. |
| | 1065 | -- sourceschema text - Database schema of the source geometry table. |
| | 1066 | -- sourcetable text - Source geometry table name. |
| | 1067 | -- sourcegeomcolumn text - Source geometry table geometry column name. |
| | 1068 | -- maxdist double precision - Maximum distance from pixel to the source |
| | 1069 | -- when the source is too far from the pixel. |
| | 1070 | -- Pixel that exceeds it will be assigned nodatavalue. |
| | 1071 | -------------------------------------------------------------------- |
| | 1072 | DROP FUNCTION IF EXISTS ST_EuclideanDistance(width integer, |
| | 1073 | height integer, |
| | 1074 | rastulx float8, |
| | 1075 | rastuly float8, |
| | 1076 | rastscalex float8, |
| | 1077 | rastscaley float8, |
| | 1078 | rastskewx float8, |
| | 1079 | rastskewy float8, |
| | 1080 | rastsrid integer, |
| | 1081 | rastndv float8, |
| | 1082 | pixeltype text, |
| | 1083 | sourceschema text, |
| | 1084 | sourcetable text, |
| | 1085 | sourcegeomcolumn text, |
| | 1086 | double precision = -1); |
| | 1087 | CREATE OR REPLACE FUNCTION ST_EuclideanDistance(width integer, |
| | 1088 | height integer, |
| | 1089 | rastulx float8, |
| | 1090 | rastuly float8, |
| | 1091 | rastscalex float8, |
| | 1092 | rastscaley float8, |
| | 1093 | rastskewx float8, |
| | 1094 | rastskewy float8, |
| | 1095 | rastsrid integer, |
| | 1096 | rastndv float8, |
| | 1097 | pixeltype text, |
| | 1098 | sourceschema text, |
| | 1099 | sourcetable text, |
| | 1100 | sourcegeomcolumn text, |
| | 1101 | double precision = -1) |
| | 1102 | RETURNS raster AS |
| | 1103 | $$ |
| | 1104 | DECLARE |
| | 1105 | maxdist ALIAS FOR $15; |
| | 1106 | x integer; |
| | 1107 | y integer; |
| | 1108 | sourcex float8; |
| | 1109 | sourcey float8; |
| | 1110 | sourcexr integer; |
| | 1111 | sourceyr integer; |
| | 1112 | sourcesrid integer; |
| | 1113 | newx float8; |
| | 1114 | newy float8; |
| | 1115 | exesql text; |
| | 1116 | newnodatavalue float8; |
| | 1117 | newinitialvalue float8; |
| | 1118 | newrast raster; |
| | 1119 | newval float8; |
| | 1120 | newpixeltype text; |
| | 1121 | -- Assuming reference raster has only one band |
| | 1122 | band integer := 1; |
| | 1123 | geom geometry; |
| | 1124 | BEGIN |
| | 1125 | -- Assuming source point table has one SRID for all geometries |
| | 1126 | EXECUTE 'SELECT DISTINCT ST_SRID(' || sourcegeomcolumn || ') FROM ' || sourceschema || '.' || sourcetable || ';' INTO sourcesrid; |
| | 1127 | |
| | 1128 | -- Create a new empty raster using provided georeference |
| | 1129 | newrast := ST_MakeEmptyRaster(width, height, rastulx, rastuly, rastscalex,rastscaley, rastskewx, rastskewy, rastsrid); |
| | 1130 | |
| | 1131 | -- Check if the new raster is empty (width = 0 OR height = 0) |
| | 1132 | IF ST_IsEmpty(newrast) THEN |
| | 1133 | RAISE NOTICE 'ST_EuclideanDistance: Raster is empty. Returning an empty raster'; |
| | 1134 | RETURN newrast; |
| | 1135 | END IF; |
| | 1136 | |
| | 1137 | -- Set the new pixeltype |
| | 1138 | newpixeltype := pixeltype; |
| | 1139 | IF newpixeltype != '1BB' AND newpixeltype != '2BUI' AND newpixeltype != '4BUI' AND newpixeltype != '8BSI' AND newpixeltype != '8BUI' AND |
| | 1140 | newpixeltype != '16BSI' AND newpixeltype != '16BUI' AND newpixeltype != '32BSI' AND newpixeltype != '32BUI' AND newpixeltype != '32BF' AND newpixeltype != '64BF' THEN |
| | 1141 | RAISE EXCEPTION 'ST_EuclideanDistance: Invalid pixeltype "%". Aborting.', newpixeltype; |
| | 1142 | END IF; |
| | 1143 | |
| | 1144 | -- Check for notada value |
| | 1145 | newnodatavalue := rastndv; |
| | 1146 | IF newnodatavalue IS NULL OR newnodatavalue < ST_MinPossibleValue(newpixeltype) OR newnodatavalue > (-ST_MinPossibleValue(newpixeltype) - 1) THEN |
| | 1147 | 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'; |
| | 1148 | newnodatavalue := ST_MinPossibleValue(newpixeltype); |
| | 1149 | END IF; |
| | 1150 | -- We set the initial value of the resulting raster to nodata value. |
| | 1151 | -- If nodatavalue is null then the raster will be initialise to ST_MinPossibleValue |
| | 1152 | newinitialvalue := newnodatavalue; |
| | 1153 | |
| | 1154 | -- Create the raster receiving all the computed distance values. Initialize it to the new initial value. |
| | 1155 | newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue); |
| | 1156 | |
| | 1157 | -- Set raster pixel value to zero if pixel is source point |
| | 1158 | exesql := 'SELECT ' || sourcegeomcolumn || ' FROM ' || sourceschema || '.' || sourcetable || ';'; |
| | 1159 | FOR geom IN EXECUTE(exesql) LOOP |
| | 1160 | sourcex := ST_X(geom); |
| | 1161 | sourcey := ST_Y(geom); |
| | 1162 | sourcexr := ST_World2RasterCoordX(newrast, sourcex, sourcey); |
| | 1163 | sourceyr := ST_World2RasterCoordY(newrast, sourcex, sourcey); |
| | 1164 | -- If source point within raster range, then set pixel value to zero |
| | 1165 | IF sourcexr <= width AND sourceyr <= height THEN |
| | 1166 | newrast := ST_SetValue(newrast, band, sourcexr, sourceyr, 0); |
| | 1167 | END IF; |
| | 1168 | END LOOP; |
| | 1169 | |
| | 1170 | -- Scan pixels in the raster set pixel value to the computed Eculidean distance |
| | 1171 | -- to its nearest source point. |
| | 1172 | -- There is place for optimization here for a more efficient scanning method |
| | 1173 | FOR x IN 1..width LOOP |
| | 1174 | FOR y IN 1..height LOOP |
| | 1175 | newx := ST_Raster2WorldCoordX(newrast, x, y); |
| | 1176 | newy := ST_Raster2WorldCoordY(newrast, x, y); |
| | 1177 | -- If pixel is not the source point then compute Eculidean distance |
| | 1178 | IF ST_Value(newrast, band, x, y) != 0 OR ST_Value(newrast, band, x, y) IS NULL THEN |
| | 1179 | exesql := 'SELECT ST_Distance(ST_SetSRID(ST_MakePoint(' || newx || ',' || newy || '),' || sourcesrid || '), (SELECT ' || sourcegeomcolumn || ' FROM ' || sourceschema || '.' || sourcetable || ' ORDER BY ' || sourcegeomcolumn || ' <-> ST_SetSRID(ST_MakePoint(' || newx || ',' || newy || '),' || sourcesrid || ') LIMIT 1));'; |
| | 1180 | EXECUTE exesql INTO newval; |
| | 1181 | -- If max distance specified |
| | 1182 | IF maxdist != -1 THEN |
| | 1183 | -- Any distance exceeds max distance is assigned nodata value |
| | 1184 | IF newval > maxdist THEN |
| | 1185 | newval = newnodatavalue; |
| | 1186 | END IF; |
| | 1187 | END IF; |
| | 1188 | END IF; |
| | 1189 | newrast := ST_SetValue(newrast, band, x, y, newval); |
| | 1190 | END LOOP; |
| | 1191 | END LOOP; |
| | 1192 | |
| | 1193 | RETURN newrast; |
| | 1194 | END; |
| | 1195 | $$ |
| | 1196 | LANGUAGE 'plpgsql'; |
| | 1197 | |
| | 1198 | -- test |
| | 1199 | -- CREATE TABLE test_eudist_2 AS (SELECT 1 AS rid,ST_EuclideanDistance(10, 10, 0, 0, 9, 9, 0, 0, 4326,-1,'32BF','public','test_geom','the_geom') AS rast FROM test_raster); |
| | 1200 | -- CREATE TABLE test_eudist_22 AS (SELECT 1 AS rid,ST_EuclideanDistance(10, 10, 0, 0, 9, 9, 0, 0, 4326,-1,'32BF','public','test_geom','the_geom',4954880) AS rast FROM test_raster); |
| | 1201 | |
| | 1202 | |
| | 1203 | -------------------------------------------------------------------- |
| | 1204 | -- ST_EuclideanDistance - (for one raster) Return a raster generated from a one band reference raster |
| | 1205 | -- which values are the Euclidean distance |
| | 1206 | -- from pixel to the points (for now only concerning point geometries) |
| | 1207 | -- from the input source table of geometries |
| | 1208 | -- (points for now, lines or polygons in the future). |
| | 1209 | -- Arguments |
| | 1210 | -- refrast raster - Reference raster align with which Euclidean distance raster will be created. |
| | 1211 | -- pixeltype text - Pixeltype assigned to the resulting raster. By default |
| | 1212 | -- to be the pixeltype of the reference raster if specified. |
| | 1213 | -- sourceschema text - Database schema of the source geometry table. |
| | 1214 | -- sourcetable text - Source geometry table name. |
| | 1215 | -- sourcegeomcolumn text - Source geometry table geometry column name. |
| | 1216 | -- maxdist double precision - Maximum distance from pixel to the source |
| | 1217 | -- when the source is too far from the pixel. |
| | 1218 | -- Pixel that exceeds it will be assigned nodatavalue. |
| | 1219 | -------------------------------------------------------------------- |
| | 1220 | DROP FUNCTION IF EXISTS ST_EuclideanDistance(refrast raster, |
| | 1221 | pixeltype text, |
| | 1222 | sourceschema text, |
| | 1223 | sourcetable text, |
| | 1224 | sourcegeomcolumn text, |
| | 1225 | double precision = -1); |
| | 1226 | CREATE OR REPLACE FUNCTION ST_EuclideanDistance(refrast raster, |
| | 1227 | pixeltype text, |
| | 1228 | sourceschema text, |
| | 1229 | sourcetable text, |
| | 1230 | sourcegeomcolumn text, |
| | 1231 | double precision = -1) |
| | 1232 | RETURNS raster AS |
| | 1233 | $$ |
| | 1234 | DECLARE |
| | 1235 | maxdist ALIAS FOR $6; |
| | 1236 | width integer; |
| | 1237 | height integer; |
| | 1238 | x integer; |
| | 1239 | y integer; |
| | 1240 | sourcex float8; |
| | 1241 | sourcey float8; |
| | 1242 | sourcexr integer; |
| | 1243 | sourceyr integer; |
| | 1244 | sourcesrid integer; |
| | 1245 | newsrid integer; |
| | 1246 | newx float8; |
| | 1247 | newy float8; |
| | 1248 | exesql text; |
| | 1249 | newnodatavalue float8; |
| | 1250 | newinitialvalue float8; |
| | 1251 | newrast raster; |
| | 1252 | newval float8; |
| | 1253 | newpixeltype text; |
| | 1254 | -- Assuming reference raster has only one band |
| | 1255 | band integer := 1; |
| | 1256 | geom geometry; |
| | 1257 | BEGIN |
| | 1258 | -- Check if reference raster is NULL |
| | 1259 | IF refrast IS NULL THEN |
| | 1260 | RAISE NOTICE 'ST_EuclideanDistance: Reference raster is NULL. Returning NULL'; |
| | 1261 | RETURN NULL; |
| | 1262 | END IF; |
| | 1263 | |
| | 1264 | width := ST_Width(refrast); |
| | 1265 | height := ST_Height(refrast); |
| | 1266 | newsrid := ST_SRID(refrast); |
| | 1267 | -- Assuming source point table has one SRID for all geometries |
| | 1268 | EXECUTE 'SELECT DISTINCT ST_SRID(' || sourcegeomcolumn || ') FROM ' || sourceschema || '.' || sourcetable || ';' INTO sourcesrid; |
| | 1269 | |
| | 1270 | -- Create a new empty raster having the same georeference as the provided raster |
| | 1271 | newrast := ST_MakeEmptyRaster(width, height, ST_UpperLeftX(refrast), ST_UpperLeftY(refrast), ST_ScaleX(refrast), ST_ScaleY(refrast), ST_SkewX(refrast), ST_SkewY(refrast), newsrid); |
| | 1272 | |
| | 1273 | -- Check if the new raster is empty (width = 0 OR height = 0) |
| | 1274 | IF ST_IsEmpty(newrast) THEN |
| | 1275 | RAISE NOTICE 'ST_EuclideanDistance: Raster is empty. Returning an empty raster'; |
| | 1276 | RETURN newrast; |
| | 1277 | END IF; |
| | 1278 | |
| | 1279 | -- Set the new pixeltype |
| | 1280 | newpixeltype := pixeltype; |
| | 1281 | -- If pixeltype not specified then use the pixeltype of the reference raster |
| | 1282 | IF newpixeltype IS NULL THEN |
| | 1283 | newpixeltype := ST_BandPixelType(refrast, band); |
| | 1284 | ELSIF newpixeltype != '1BB' AND newpixeltype != '2BUI' AND newpixeltype != '4BUI' AND newpixeltype != '8BSI' AND newpixeltype != '8BUI' AND |
| | 1285 | newpixeltype != '16BSI' AND newpixeltype != '16BUI' AND newpixeltype != '32BSI' AND newpixeltype != '32BUI' AND newpixeltype != '32BF' AND newpixeltype != '64BF' THEN |
| | 1286 | RAISE EXCEPTION 'ST_EuclideanDistance: Invalid pixeltype "%". Aborting.', newpixeltype; |
| | 1287 | END IF; |
| | 1288 | |
| | 1289 | -- Check for notada value |
| | 1290 | newnodatavalue := ST_BandNodataValue(refrast, band); |
| | 1291 | IF newnodatavalue IS NULL OR newnodatavalue < ST_MinPossibleValue(newpixeltype) OR newnodatavalue > (-ST_MinPossibleValue(newpixeltype) - 1) THEN |
| | 1292 | 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'; |
| | 1293 | newnodatavalue := ST_MinPossibleValue(newpixeltype); |
| | 1294 | END IF; |
| | 1295 | -- We set the initial value of the resulting raster to nodata value. |
| | 1296 | -- If nodatavalue is null then the raster will be initialise to ST_MinPossibleValue |
| | 1297 | newinitialvalue := newnodatavalue; |
| | 1298 | |
| | 1299 | -- Create the raster receiving all the computed distance values. Initialize it to the new initial value. |
| | 1300 | newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue); |
| | 1301 | |
| | 1302 | -- Set raster pixel value to zero if pixel is source point |
| | 1303 | exesql := 'SELECT ' || sourcegeomcolumn || ' FROM ' || sourceschema || '.' || sourcetable || ';'; |
| | 1304 | FOR geom IN EXECUTE(exesql) LOOP |
| | 1305 | sourcex := ST_X(geom); |
| | 1306 | sourcey := ST_Y(geom); |
| | 1307 | sourcexr := ST_World2RasterCoordX(newrast, sourcex, sourcey); |
| | 1308 | sourceyr := ST_World2RasterCoordY(newrast, sourcex, sourcey); |
| | 1309 | -- If source point within raster range, then set pixel value to zero |
| | 1310 | IF sourcexr <= width AND sourceyr <= height THEN |
| | 1311 | newrast := ST_SetValue(newrast, band, sourcexr, sourceyr, 0); |
| | 1312 | END IF; |
| | 1313 | END LOOP; |
| | 1314 | -- Scan pixels in the raster set pixel value to the computed Eculidean distance |
| | 1315 | -- to its nearest source point. |
| | 1316 | -- There is place for optimization here for a more efficient scanning method |
| | 1317 | FOR x IN 1..width LOOP |
| | 1318 | FOR y IN 1..height LOOP |
| | 1319 | newx := ST_Raster2WorldCoordX(newrast, x, y); |
| | 1320 | newy := ST_Raster2WorldCoordY(newrast, x, y); |
| | 1321 | -- If pixel is not the source point then compute Eculidean distance |
| | 1322 | IF ST_Value(newrast, band, x, y) != 0 OR ST_Value(newrast, band, x, y) IS NULL THEN |
| | 1323 | exesql := 'SELECT ST_Distance(ST_SetSRID(ST_MakePoint(' || newx || ',' || newy || '),' || sourcesrid || '), (SELECT ' || sourcegeomcolumn || ' FROM ' || sourceschema || '.' || sourcetable || ' ORDER BY ' || sourcegeomcolumn || ' <-> ST_SetSRID(ST_MakePoint(' || newx || ',' || newy || '),' || sourcesrid || ') LIMIT 1));'; |
| | 1324 | EXECUTE exesql INTO newval; |
| | 1325 | -- If max distance specified |
| | 1326 | IF maxdist != -1 THEN |
| | 1327 | -- Any distance exceeds max distance is assigned nodata value |
| | 1328 | IF newval > maxdist THEN |
| | 1329 | newval := newnodatavalue; |
| | 1330 | END IF; |
| | 1331 | END IF; |
| | 1332 | END IF; |
| | 1333 | newrast := ST_SetValue(newrast, band, x, y, newval); |
| | 1334 | END LOOP; |
| | 1335 | END LOOP; |
| | 1336 | |
| | 1337 | RETURN newrast; |
| | 1338 | END; |
| | 1339 | $$ |
| | 1340 | LANGUAGE 'plpgsql'; |
| | 1341 | |
| | 1342 | -- test |
| | 1343 | -- CREATE TABLE test_eudist AS (SELECT 1 AS rid,ST_EuclideanDistance(rast,NULL,'public','test_geom','the_geom') AS rast FROM test_raster); |
| | 1344 | -- CREATE TABLE test_eudist2 AS (SELECT 1 AS rid,ST_EuclideanDistance(rast,NULL,'public','test_geom','the_geom',4954880) AS rast FROM test_raster); |
| | 1345 | }}} |
