Sunday, 26 June 2016

grass - How to identify line intersection in QGIS when I have more than 2 lines?


I am trying to identify where overlaps exist between multiple shapefiles. I would use the intersect tool, but it only allows me to intersect two layers.


I have over one hundred individual line shapefiles and I want line segments that overlap.



Any ideas?


I'm using QGIS 2.0.1 and GRASS for the creation of the shapefiles.



Answer



The idea is to merge all your shapefiles containing linestrings. You have to handle with unique identifiers that are characterizing your linestrings definitely, according to the distinct shapefiles.


QGIS approach:




  1. Add the unique IDs to the attribute table of each shapefile if required: Open the attribute table, toggle editing and add a new column. Fill the column with a unique ID, e.g. with the field calculator.





  2. Merge all your shapefiles. See MMqgis (mmqgis > combine > merge layers) and further information here.




  3. Use the line intersection tool (vector > analysis tool > line intersection). You have to intersect the merged shapefile with itself. The result is a point shapefile (layer). In this shapefile you see all the points that are representing an intersection and both IDs which are representing the intersected linestrings. You can then sort the intersections by the IDs.




UPDATE


I have overseen that you want to have the overlapping segments. Sure, you can do step 3 as well with Vector > Geoprocessing Tools > Intersect. After that use the filter as mentioned below. Maybe the points which are symbolizing the intersection between lines are interesting for you as well.


You have to use GRASS GIS in addition, because I think it is impossible in QGIS to add columns with unique IDs to hundreds of shapefiles.


I tried this approach with a small sample of data.



This is the attribute (point) table from the intersection (Analysis Tool) of three linestring layers.


enter image description here


On this picture you see which line intersects which line.


enter image description here


The one intersection (blue, blue) has not a resulting point, because it is the same linestring.


What you finally need:


Here is the picture of the overlaps after the intersection with Geoprocessing Tool. I used some other IDs here.


enter image description here


You can then filter your results with right clicking on the intersection-segment-layer and clicking on filter (filter expression in the lowermost box). For instance the intersection between linestrings with id=1 and id_2=11 (in my sample data): "id" = 1 AND "id_2" = 11.


GRASS GIS approach (console based)



Open your console or terminal, start GRASS and open your mapset




  1. Import your shapefiles into your GRASS mapset. Before, create a directory for only(!) all your linestring layers.


    for i in *.shp; do v.in.ogr -o dsn=/your_shapefile_directory layer=${i%.*} output=${i%.*}; done


  2. Adding a new column to your linestring layers in GRASS if required. The parameter pattern stands for the name of your linestring layers. In my case I have the linestring layers line1, line2 and line3. With the expression "line*" you select them all. Use pattern="*" to select all layers without any taxonomy.


    for i in `g.mlist type=vect pattern="line*"`; do v.db.addcol map=$i columns="gid int"; done


    The new column is gid from type integer.




  3. Now you have to fill the columns (gid) for every layer


    val=1
    for i in `g.mlist type=vect pattern="line*"`; do v.db.update $i column=gid value=$val; val=$((val+1)); done

    First, define the val(hit ENTER), then run the loop





  4. At least, you merge your layers


    names=""
    for i in `g.mlist type=vect pattern="line*"`; do names="$i,$names"; done
    v.patch -e input=${names%,} output=merged_lines

    First you define names, which will be the list for the input. Then you run the loop. The resulting list is line1,line2,line3,. After that you merge the list avoiding the last comma. You have to specify -e for merging the attribute table as well.




  5. Load the merged shapefile from grass into QGIS and do the final updated step 3 from the QGIS approach (at the beginning of this answer).





No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...