DROP TYPE IF EXISTS otm_al_grid; CREATE TYPE otm_al_grid AS (x DOUBLE PRECISION, y DOUBLE PRECISION, is_in_area BOOLEAN, is_near_axis BOOLEAN, dist_to_border DOUBLE PRECISION, dist_to_middle DOUBLE PRECISION, way_to_middle INTEGER ); CREATE OR REPLACE FUNCTION arealabel(myosm_id IN BIGINT,myway IN GEOMETRY) RETURNS GEOMETRY AS $$ DECLARE bbox GEOMETRY; tmpway GEOMETRY; linestring1 VARCHAR; linestring2 VARCHAR; t VARCHAR; x DOUBLE PRECISION; y DOUBLE PRECISION; xmin DOUBLE PRECISION; xmax DOUBLE PRECISION; ymin DOUBLE PRECISION; ymax DOUBLE PRECISION; i INTEGER; j INTEGER; i1 INTEGER; j1 INTEGER; imin INTEGER; jmin INTEGER; imax INTEGER; jmax INTEGER; p INTEGER; k DOUBLE PRECISION; f DOUBLE PRECISION; g DOUBLE PRECISION; m DOUBLE PRECISION; changes INTEGER; gridsize CONSTANT INTEGER:=30; grid otm_al_grid[900]; gridpoint otm_al_grid; targetpoint otm_al_grid; gridwidth DOUBLE PRECISION; lastindex INTEGER; middleindex INTEGER; endindex1 INTEGER; endindex2 INTEGER; BEGIN -- -- get bounding box, adjust grid (longer side is 'gridsize') to same grid with for x and y -- bbox:=ST_SetSRID(ST_Envelope(myway),3857); xmin:=ST_XMin(bbox);xmax:=ST_XMax(bbox); ymin:=ST_YMin(bbox);ymax:=ST_YMax(bbox); x=xmax-xmin;y:=ymax-ymin; IF (x>y) THEN gridwidth:=x/(gridsize-1); ELSE gridwidth:=y/(gridsize-1); END IF; x:=xmin;i:=0;j:=0; -- -- initialize grid point with coordinates, check which grid point is in the area -- WHILE (x<=xmax) LOOP y:=ymin;j:=0; WHILE (y<=ymax) LOOP gridpoint.x=x; gridpoint.y=y; gridpoint.is_in_area=ST_Within(st_Expand(ST_SetSRID(ST_Point(x,y),3857),gridwidth/2),ST_SetSRID(myway,3857)); gridpoint.is_near_axis:=false; gridpoint.dist_to_border:=0; gridpoint.dist_to_middle:=0; gridpoint.way_to_middle:=0; grid[i*gridsize+j]:=gridpoint; y:=y+gridwidth; j:=j+1; END LOOP; x:=x+gridwidth; i:=i+1; END LOOP; -- -- calc distance to borders -- FOR i IN 1..(gridsize-2) LOOP FOR j IN 1..(gridsize-2) LOOP gridpoint=grid[i*gridsize+j]; IF (gridpoint.is_in_area) THEN gridpoint.dist_to_border=1+LEAST(grid[(i-1)*gridsize+j].dist_to_border,grid[(i)*gridsize+(j-1)].dist_to_border); END IF; grid[i*gridsize+j]:=gridpoint; END LOOP; END LOOP; FOR i IN REVERSE (gridsize-2)..1 LOOP FOR j IN REVERSE (gridsize-2)..1 LOOP gridpoint=grid[i*gridsize+j]; IF (gridpoint.is_in_area) THEN gridpoint.dist_to_border=LEAST(gridpoint.dist_to_border,1+LEAST(grid[(i+1)*gridsize+j].dist_to_border,grid[(i)*gridsize+(j+1)].dist_to_border)); END IF; grid[i*gridsize+j]:=gridpoint; END LOOP; END LOOP; -- -- mark points with >=3 neighbours, we don't need these points near the border -- middleindex:=0;p:=0;imin=gridsize;jmin=gridsize;imax:=-1;jmax:=-1; FOR i IN 1..(gridsize-2) LOOP FOR j IN 1..(gridsize-2) LOOP gridpoint=grid[i*gridsize+j]; IF (gridpoint.is_in_area) THEN i1:=0; IF (gridpoint.dist_to_border>=grid[(i-1)*gridsize+j].dist_to_border) THEN i1:=i1+1; END IF; IF (gridpoint.dist_to_border>=grid[(i+1)*gridsize+j].dist_to_border) THEN i1:=i1+1; END IF; IF (gridpoint.dist_to_border>=grid[(i)*gridsize+j-1].dist_to_border) THEN i1:=i1+1; END IF; IF (gridpoint.dist_to_border>=grid[(i)*gridsize+j+1].dist_to_border) THEN i1:=i1+1; END IF; IF (i1>2) THEN gridpoint.is_near_axis:=true; if(imin>i) THEN imin:=i; END IF; if(jmin>j) THEN jmin:=j; END IF; if(imaxp) THEN p:=gridpoint.dist_to_border;middleindex:=(i)*gridsize+j; END IF; END IF; END IF; END LOOP; END LOOP; -- -- calc distance to the center point =(one of) the point(s) with the max dist to the border -- mark (one of the) node(s) with max. distance -- gridpoint=grid[middleindex]; gridpoint.dist_to_middle:=1; grid[middleindex]:=gridpoint; changes:=1;g:=0.0;m:=0.0; WHILE (changes>0) LOOP changes:=0;m:=m+1; FOR i IN imin..imax LOOP FOR j IN jmin..jmax LOOP gridpoint=grid[i*gridsize+j]; IF ((gridpoint.is_near_axis) AND (gridpoint.dist_to_middle>0.0)) THEN FOR i1 IN (i-1)..(i+1) LOOP FOR j1 IN (j-1)..(j+1) LOOP IF ((i1!=i)OR(j1!=j)) THEN targetpoint:=grid[i1*gridsize+j1]; IF (targetpoint.is_near_axis) THEN IF ((i=i1) OR (j=j1)) THEN f:=1.0; ELSE f:=1.4142 ; END IF; IF ((targetpoint.dist_to_middle=0.0) OR (targetpoint.dist_to_middle>gridpoint.dist_to_middle+f)) THEN targetpoint.dist_to_middle:=gridpoint.dist_to_middle+f; targetpoint.way_to_middle:=i*gridsize+j; grid[i1*gridsize+j1]:=targetpoint; if(targetpoint.dist_to_middle>g) THEN g:=targetpoint.dist_to_middle;endindex1:=i1*gridsize+j1; END IF; changes:=changes+1; END IF; END IF; END IF; END LOOP; END LOOP; END IF; END LOOP; END LOOP; FOR i IN REVERSE imax..imin LOOP FOR j IN REVERSE jmax..jmin LOOP gridpoint=grid[i*gridsize+j]; IF ((gridpoint.is_near_axis) AND (gridpoint.dist_to_middle>0.0)) THEN FOR i1 IN (i-1)..(i+1) LOOP FOR j1 IN (j-1)..(j+1) LOOP IF ((i1!=i)OR(j1!=j)) THEN targetpoint:=grid[i1*gridsize+j1]; IF (targetpoint.is_near_axis) THEN IF ((i=i1) OR (j=j1)) THEN f:=1.0; ELSE f:=1.4142 ; END IF; IF ((targetpoint.dist_to_middle=0.0) OR (targetpoint.dist_to_middle>gridpoint.dist_to_middle+f)) THEN targetpoint.dist_to_middle:=gridpoint.dist_to_middle+f; targetpoint.way_to_middle:=i*gridsize+j; grid[i1*gridsize+j1]:=targetpoint; if(targetpoint.dist_to_middle>g) THEN g:=targetpoint.dist_to_middle;endindex1:=i1*gridsize+j1; END IF; changes:=changes+1; END IF; END IF; END IF; END LOOP; END LOOP; END IF; END LOOP; END LOOP; END LOOP; RAISE NOTICE ' MID to 1st node ID: % rounds: % % %',myosm_id,m,middleindex,endindex1; -- -- Build linestring 1.st node to middle -- i:=endindex1;lastindex:=i; linestring1:='LINESTRING('; WHILE (i!=middleindex) LOOP gridpoint=grid[i]; linestring1:=linestring1 || gridpoint.x || ' ' || gridpoint.y || ','; lastindex:=i; i:=gridpoint.way_to_middle; END LOOP; gridpoint=grid[i]; linestring1:=linestring1 || gridpoint.x || ' ' || gridpoint.y || ')'; -- -- Find node on the other side of the middle. To avoid sharp edges at the middle point -- search only in the opposite direction of the last part of the way. -- if((lastindex/gridsize)>(i/gridsize)) THEN imax=i/gridsize; END IF; if((lastindex/gridsize)<(i/gridsize)) THEN imin=i/gridsize; END IF; if((lastindex%gridsize)>(i%gridsize)) THEN jmax=i%gridsize; END IF; if((lastindex%gridsize)<(i%gridsize)) THEN jmin=i%gridsize; END IF; g:=0; FOR i IN imin..imax LOOP FOR j IN jmin..jmax LOOP gridpoint=grid[i*gridsize+j]; IF ((i*gridsize+j!=endindex1) AND (gridpoint.dist_to_middle>g)) THEN endindex2:=i*gridsize+j; g:=gridpoint.dist_to_middle; END IF; END LOOP; END LOOP; -- -- Build linestring 2nd node to middle -- i:=endindex2; linestring2:='LINESTRING('; WHILE (i!=middleindex) LOOP gridpoint=grid[i]; linestring2:=linestring2 || gridpoint.x || ' ' || gridpoint.y || ','; i:=gridpoint.way_to_middle; END LOOP; gridpoint=grid[i]; linestring2:=linestring2 || gridpoint.x || ' ' || gridpoint.y || ')'; -- -- DEBUG: insert nodes and lines -- FOR i IN 0..(gridsize-1) LOOP FOR j IN 0..(gridsize-1) LOOP IF (grid[(i)*gridsize+j].is_in_area) THEN gridpoint=grid[i*gridsize+j]; x:=gridpoint.x; y:=gridpoint.y; tmpway:=ST_SetSRID(ST_Point(x,y),3857); t:=grid[(i)*gridsize+j].dist_to_border::TEXT; IF (gridpoint.is_near_axis) THEN INSERT into planet_osm_point (osm_id,way,man_made,amenity,name) values (myosm_id,tmpway,'lakepoint','nearaxis',t); ELSE INSERT into planet_osm_point (osm_id,way,man_made,amenity,name) values (myosm_id,tmpway,'lakepoint','notaxis',t); END IF; END IF; END LOOP; END LOOP; IF (endindex1!=middleindex) THEN INSERT into planet_osm_line (osm_id,way,name,man_made,amenity) values (myosm_id,(ST_GeomFromText(linestring1,3857)),'XXX','lakeaxis','1'); -- INSERT into planet_osm_line (osm_id,way,name,man_made,amenity) values (myosm_id,st_simplify((ST_GeomFromText(linestring1,3857)),gridwidth*1.4142),'XXX','lakeaxis','1'); ELSE RAISE NOTICE 'unable to find way in area %',myosm_id; END IF; IF (endindex2!=middleindex) THEN INSERT into planet_osm_line (osm_id,way,name,man_made,amenity) values (myosm_id,(ST_GeomFromText(linestring2,3857)),'XXX','lakeaxis','2'); -- INSERT into planet_osm_line (osm_id,way,name,man_made,amenity) values (myosm_id,st_simplify((ST_GeomFromText(linestring2,3857)),gridwidth*1.4142),'XXX','lakeaxis','2'); ELSE RAISE NOTICE 'unable to find way in area %',myosm_id; END IF; RETURN NULL; END; $$ LANGUAGE plpgsql; delete from planet_osm_point where man_made='lakepoint'; delete from planet_osm_line where man_made='lakeaxis'; SELECT osm_id,name,arealabel(osm_id,way) from planet_osm_polygon where name is not null and "natural"='water';