I am attempting to select lines that radiate from a point using v.select in GRASS GIS 7.6 (also tested with 7.4.4) in Windows 10.

I start out with 48 lines radiating from a point every 7.5 degrees. These are then clipped by the vector representing the sea (as seen in figure 1), and I subsequently want to retrieve the lines that are still touching the point.

enter image description here

However, running the following command only selects 11 of these lines (figure 2).

v.select ainput=temp_radials_clipped atype=line binput=random_point btype=point output=lines_vselect operator=overlap

enter image description here

All the lines share the exact same origin point at the coordinates of random_point. I have tried all the different operators for v.select, and only 'overlap' returns any lines at all (this also includes 'intersects' which is to be the GEOS equivalent of 'overlap'). I have trouble understanding why v.select with 'overlap' should return any lines at all if it is not returning all of the lines touching the point.

I am also open to alternatives to v.select to achieve this.

  • 1
    It's possible that computer rounding errors are changing the clipped lines very slightly, then there is no intersection. Try buffering the random point by a small amout (maybe 0.1 m will be enough) then try to intersect the radial lines with that buffer. – Micha Jun 14 '19 at 19:55
  • Thank you, that seems to answer it. I can push the buffer as far down as 0.0000001 m and still get the correct output using v.select with the 'intersects' operator. Would best practice now be that you post your comment as an answer and I accept it as correct? Apologies, I'm quite new to this. – humperderp Jun 15 '19 at 10:41
  • 1
    No accept necessary, since I didn't really give a full answer, just a hint. But in general, the issue of intersecting/overlapping lines and points in GIS comes up often. The issues of computer rounding errors makes this very tricky. I would go so far as to say: never depend on intersection with a point. Hope that SE GIS continues to be of help. – Micha Jun 15 '19 at 14:25

This answer is largely based on the comments by Micha, posted to the original question. These are reiterated here so that the question can be marked as answered.

The problem seems to arise from computer rounding error that results in the lines being ever so slightly off the point. One work-around is creating a small buffer around the point and running v.select using the 'intersects' operator on the lines and buffer instead of on the point. Creating a buffer of 0.1 m works fine for my purposes, but it appears this can be reduced even further, should that be necessary.

Below is a more long-winded alternative via a python script that has worked so far. This is based on the first map provided in the question. Apart from the number of steps required, this seems less desirable than the 'v.select-on-buffer'-option as it is also based on the intersection of lines with a point. This might therefore conceivably run into the same problem of rounding error at some point (cf. comment by Micha).

# Create the new vector target_lines to hold the relevant lines 
grass.run_command('v.edit', map = 'target_lines', tool = 'create')  

# Retrieve the coordinates of the point
report = grass.read_command('v.report', map = 'random_point', option= 'coor')
ptCoords = report.split('\n')[1].split('|')[1:3]

# Read the IDs of the features of temp_radials_clipped that are at the coordinates of the point
lineIds = grass.read_command('v.edit', map = 'temp_radials_clipped', tool = 'select', 
                             type = 'line', coords = ptCoords)

# Copy these features over to target_lines
grass.run_command('v.edit', map = 'target_lines', tool = 'copy', bgmap = 'temp_radials_clipped',
                  ids = lineIds)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.