Grid snapping only in SQL!

run.sql

CREATE OR REPLACE FUNCTION axh_miles_test(input GEOMETRY)
  RETURNS GEOMETRY AS $$
DECLARE
  -- sql TEXT;
  first GEOMETRY;
  next GEOMETRY;
  back GEOMETRY;
  last GEOMETRY;
  invert GEOMETRY;
  linestring GEOMETRY;
  backstring GEOMETRY;
  -- mem GEOMETRY;
  nearest GEOMETRY;
  closest GEOMETRY;
  -- start BOOLEAN := TRUE;
  ret axh_miles_test_ret%rowtype;
  ret_id INT[];
  ret_geom GEOMETRY[];
  ret_line GEOMETRY[];
  size INT;
  curr INT;
  top GEOMETRY[];
  bot GEOMETRY[];
  -- allret axh_miles_test_ret%rowtype[];
  --CREATE TYPE axh_miles_test_ret AS (i INT, id INT, geom GEOMETRY, line GEOMETRY);
BEGIN

   first := ST_PointN(input, 1);
   last := ST_PointN(ST_Reverse(input), 1);


   SELECT INTO nearest line_interpolate_point(ST_GeometryN(the_geom,1), line_locate_point(ST_GeometryN(the_geom,1), first)) FROM new_york_osm_line_clean WHERE z_order > 2 AND the_geom IS NOT NULL ORDER BY the_geom <-> first LIMIT 1;

   -- top := top || first;
   select into top array_agg(i) from (values (first), (nearest)) n(i);
   -- top := top || nearest;
   
   -- init our linestring, we'll skip the first array point later
   linestring := ST_MakeLine(first, nearest);

   -- DO IT BACKWARDS 
   invert = ST_Reverse(input);

   SELECT INTO closest line_interpolate_point(ST_GeometryN(the_geom,1), line_locate_point(ST_GeometryN(the_geom,1), last)) FROM new_york_osm_line_clean WHERE z_order > 2 AND the_geom IS NOT NULL ORDER BY the_geom <-> last LIMIT 1;

   select into bot array_agg(i) from (values (last), (closest)) n(i);

   backstring := ST_MakeLine(last, closest);

    size := 200;
    curr := 1;
    WHILE curr <= size LOOP

      -- Do it from the top
      -- first := top[1:array_upper(top, 1) - 1];
      first := ST_Line_Interpolate_Point(ST_MakeLine(nearest, last), (curr+0.0)/(2.0*size+1.0));
      IF curr <(size/4) THEN 
        SELECT INTO nearest line_interpolate_point(ST_GeometryN(the_geom,1), line_locate_point(ST_GeometryN(the_geom,1), first)) FROM new_york_osm_line_clean WHERE z_order > 2 AND the_geom IS NOT NULL ORDER BY the_geom <-> first LIMIT 1;
      ELSE
        SELECT INTO nearest line_interpolate_point(ST_GeometryN(the_geom,1), line_locate_point(ST_GeometryN(the_geom,1), first)) FROM new_york_osm_line_clean WHERE z_order > 3 AND the_geom IS NOT NULL ORDER BY the_geom <-> first LIMIT 1;
      END IF;
      -- linestring = ST_AddPoint(linestring, nearest);
      -- top := array_append(top, nearest);
       linestring = ST_AddPoint(linestring, nearest);  
      curr := curr + 1;
    END LOOP;

    curr := 1;
    WHILE curr <= size LOOP
      -- And reverse from the bottom
      last = ST_Line_Interpolate_Point(ST_MakeLine(closest, first), (curr+0.0)/(2.0*size+1.0));
      IF curr <(size/4) THEN 
        SELECT INTO closest line_interpolate_point(ST_GeometryN(the_geom,1), line_locate_point(ST_GeometryN(the_geom,1), last)) FROM new_york_osm_line_clean WHERE z_order > 2 AND the_geom IS NOT NULL ORDER BY the_geom <-> last LIMIT 1;
      ELSE
        SELECT INTO closest line_interpolate_point(ST_GeometryN(the_geom,1), line_locate_point(ST_GeometryN(the_geom,1), last)) FROM new_york_osm_line_clean WHERE z_order > 3 AND the_geom IS NOT NULL ORDER BY the_geom <-> last LIMIT 1;
      END IF;

       backstring = ST_AddPoint(backstring, closest);  
      bot := array_append(bot, closest);
      curr := curr + 1;
    END LOOP;

    -- linestring = ST_LineMerge(ST_Collect(linestring, ST_Reverse(backstring)));
    curr := array_upper(bot, 1);
    WHILE curr > 0 LOOP
       linestring = ST_AddPoint(linestring, bot[curr]);  
       curr := curr - 1;
    END LOOP;
    -- linestring = ST_AddPoint(linestring, last);

   RETURN linestring;
END
$$ language plpgsql;