[GIS] Snap to grid with sf in R

rsfsnapping

Is there a snap to grid function available in sf for R similar to ST_SnapToGrid(geometry geomA, float size)in PostGIS? It is mentioned as a solution to non-noded intersection problem for PostGIS and am looking for a solution in R. Using other snap options with just the polygon lead to geometry errors (even with small tolerance).

Best Answer

Rounding coordinates may be equivalent to snapping all of the shapes coordinates to a regular grid.

Before -

> pnt = st_point(c(0,0))
> pol = st_buffer(pnt, 1)
> plot(pol)

After -

> pol[[1]] = round(pol[[1]], 1)
> plot(pol)

Related Solutions

[GIS] Create raster with grid cells values based on area covered by polygons

Since you put GRASS into the tags, here's a solution based on GRASS:

First, you need to know exactly what coordinate system the original data is in (as always with GRASS). I see that the *.prj file contains "TRANSVERSE MERCATOR" but it's not one of the standard UTM zones. Since you have not mentioned the EPSG code of this data, here's the proj4 string for creating a new matching GRASS Location:

+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=WGS84 +units=m +no_defs

Now start GRASS, create a LOCATION (based on the above) and a MAPSET.

And now import the shapefile into GRASS by running:

v.in.ogr input="sp_data.shp" output=sp_data

~~Add a column for the area of each feature, and calculate its area:~~

~~v.db.addcolumn map=sp_data columns="area_sqm INT" v.to.db sp_data option=area column=area_sqm unit=meters~~

~~Now convert to raster, using the area_sqm column as the raster value.~~ But first set the GRASS computational region:

g.region -p vector=sp_data res=100

Now create a vector grid of polygons at the raster resolution:

v.mkgrid map=grid

~~Intersect with the forest polygons, and calculate areas of these intersect polygons:~~ Intersect with the forest polygons, dissolve the polygons inside the same grid cell, calculate areas of these intersect polygons, and move the area value to the grid cells:

v.overlay ain=grid bin=sp_data operator=and output=forest_grid_intersect
v.dissolve --overwrite input=forest_grid_intersect@jn column=a_cat output=forest_grid_intersect_dissolve
v.db.addtable map=forest_grid_intersect_dissolve columns="area_sqm INT"
v.to.db forest_grid_intersect_dissolve option=area column=area_sqm unit=meters 
v.db.join map=grid@jn column=cat other_table=forest_grid_intersect_dissolve other_column=cat

~~v.db.addcolumn map=forest_grid_intersect columns="area_sqm INT" v.to.db forest_grid_intersect option=area column=area_sqm unit=meters~~

and convert to raster, using the area attribute as the raster value:

v.to.rast input=grid type=area output=sp_data_rast use=attr attribute_column=area_sqm

And you're done. Hope I have it right this time :-) .

QGIS – How to Snap a Road Network to a Hexagonal Grid in QGIS

My solution involves a PyQGIS script that is faster and more effective than a workflow involving snapping (I gave it a try too). Using my algorithm I've obtained these results:

You can run the following code snippets in sequence from within QGIS (in the QGIS Python console). At the end you get a memory layer with the snapped routes loaded into QGIS.

The only prerequisite is to create a multipart road Shapefile (use Processing->Singleparts to multipart, I used the field fictitiuos as Unique ID field parameter). This will give us a roads_multipart.shp file with a single feature.

Here is the algorithm explained:

Get the nearest hexagon sides where routes cross. For each hexagon we create 6 triangles between each pair of neighbour vertices and the corresponding centroid. If any road intersects a triangle, the segment shared by the hexagon and the triangle is added to the final snapped route. This is the heavier part of the whole algorithm, it takes 35 seconds running on my machine. In the first two lines there are 2 Shapefile paths, you should adjust them to fit your own file paths.

hexgrid = QgsVectorLayer("/docs/borrar/hex_grid_question/layers/normal-hexgrid.shp", "hexgrid", "ogr")
roads = QgsVectorLayer("/docs/borrar/hex_grid_question/layers/roads_multipart.shp", "roads", "ogr")  # Must be multipart!

roadFeat = roads.getFeatures().next() # We just have 1 geometry
road = roadFeat.geometry() 
indicesHexSides = ((0,1), (1,2), (2,3), (3,4), (4,5), (5,0))

epsilon = 0.01
# Function to compare whether 2 segments are equal (even if inverted)
def isSegmentAlreadySaved(v1, v2):
    for segment in listSegments:        
        p1 = QgsPoint(segment[0][0], segment[0][1])
        p2 = QgsPoint(segment[1][0], segment[1][1])
        if v1.compare(p1, epsilon) and v2.compare(p2, epsilon) \
            or v1.compare(p2, epsilon) and v2.compare(p1, epsilon):
            return True
    return False

# Let's find the nearest sides of hexagons where routes cross
listSegments = []
for hexFeat in hexgrid.getFeatures():
    hex = hexFeat.geometry()
    if hex.intersects( road ):
        for side in indicesHexSides:
            triangle = QgsGeometry.fromPolyline([hex.centroid().asPoint(), hex.vertexAt(side[0]), hex.vertexAt(side[1])])
            if triangle.intersects( road ):
                # Only append new lines, we don't want duplicates!!!
                if not isSegmentAlreadySaved(hex.vertexAt(side[0]), hex.vertexAt(side[1])): 
                    listSegments.append( [[hex.vertexAt(side[0]).x(), hex.vertexAt(side[0]).y()], [hex.vertexAt(side[1]).x(),hex.vertexAt(side[1]).y()]] )

Get rid of disconnected (or 'open') segments by using Python lists, tuples, and dictionaries. At this point, there are some disconnected segments left, i.e., segments that have one vertex disconnected but the other one connected to at least other 2 segments (see red segments in the next figure). We need to get rid of them.

# Let's remove disconnected/open segments
lstVertices = [tuple(point) for segment in listSegments for point in segment]
dictConnectionsPerVertex = dict((tuple(x),lstVertices.count(x)-1) for x in set(lstVertices))

# A vertex is not connected and the other one is connected to 2 segments
def segmentIsOpen(segment):
    return dictConnectionsPerVertex[tuple(segment[0])] == 0 and dictConnectionsPerVertex[tuple(segment[1])] >= 2 \
        or dictConnectionsPerVertex[tuple(segment[1])] == 0 and dictConnectionsPerVertex[tuple(segment[0])] >= 2

# Remove open segments
segmentsToDelete = [segment for segment in listSegments if segmentIsOpen(segment)]        
for toBeDeleted in segmentsToDelete:
    listSegments.remove( toBeDeleted )

Now we can create a vector layer from the list of coordinates and load it to the QGIS map:

# Create a memory layer and load it to QGIS map canvas
vl = QgsVectorLayer("LineString", "Snapped Routes", "memory")
pr = vl.dataProvider()
features = []
for segment in listSegments:
    fet = QgsFeature()
    fet.setGeometry( QgsGeometry.fromPolyline( [QgsPoint(segment[0][0], segment[0][1]), QgsPoint(segment[1][0], segment[1][1])] ) )
    features.append(fet)

pr.addFeatures( features )
vl.updateExtents()
QgsMapLayerRegistry.instance().addMapLayer(vl)

Another part of the result:

Should you need attributes in the snapped routes, we could use a Spatial Index to rapidly evaluate intersections (such as in https://gis.stackexchange.com/a/130440/4972 ), but that's another story.

Hope this helps!

Best Answer

Related Solutions

[GIS] Create raster with grid cells values based on area covered by polygons

QGIS – How to Snap a Road Network to a Hexagonal Grid in QGIS

Related Question