13

I'm a newbie, apologies if this is obvious/has already been asked and answered but I couldn't find anything.

I have two shapefiles: 1. an administrative boundary layer for a county in the UK known as an LSOA boundary that has 500 little zones in it 2. a flood zone.

Ideally I want to find out which of the little LSOA zones are ≥50% within the flood zone and end up with a yes/no or a 1/0 for each of the 500 LSOA zones.

But I don't know how to do this. I figured I could Join the two shapefiles, but there's no common attribute between them. Then I thought I could use the Join Attribute by Location function, which worked and shows me which LSOA are in the flood zone, but that's nearly all of them (see image 2).

I think this is an SQL problem but I don't know. I'm new to QGIS and have never used PostgreSQL.

Any help would be greatly appreciated. I can provide whatever info you lovely people need to help me out.

enter image description here

enter image description here

| improve this question | |
12

This is a relatively simple task using the geoprocessing tools included in QGIS.

  1. Calculate the area of your LSOA zones.

    • Open the LSOA layer attribute table.
    • Enable editing mode.
    • Open the field calculator.
    • Create a new field of type "Decimal number (real)" with the expression "$area".
    • Disable editing mode (saving edits).

  2. Merge the flood zone layer into a single multi-part feature.

    • Vector > Geometry Tools > Singleparts to Multipart.
    • Select "--- Merge all ---" for the Unique ID field.

  3. Intersect the LSOA zone layer with the multipart flood zone layer.

    • Vector > Geoprocessing Tools > Intersect.
    • Input layer is the LSOA zones, intersect layer is the flood zones.

  4. The resulting layer will be the parts of the LSOA zones (with the attributes from the LSOA zones layer) which overlapped with the flood zones layer. To calculate the proportion of each LSOA zone within a flood zone:

    • Calculate the area of the intersected features (as in step #1), then
    • Add another field, dividing the original (total) area by the intersected area. The result is a decimal between 0 and 1. Multiply by 100 to give a percentage.

  5. Join the original LSOA layer to the intersected layer, using the unique ID shared by both layers.

  6. Export the joined layer as a new shapefile.

  7. Delete the duplicated attributes.

Et voilà!

Without step #2, an individual feature would be created for each different flood zone feature for each LSOA feature. This probably isn't what you want if you're only interested in the total coverage for each LSOA zone. If you want to differentiate between fluvial / tidal / pluvial flooding (and the flood zone data supports it), you could convert singlepart s to multipart specifying the "TYPE" field as the Unique ID field.

| improve this answer | |
  • Thank you for your help! Much appreciated. However, I'm having some trouble. I've followed the steps. Step 3, the Intersect, took 10 hours to complete and when it was done all I got was an empty shapefile: i.imgur.com/QIM6Gtg.png Is there something I've missed? I tried to complete the process and do step 4 but there's no data to calculate the intersect area. – KJGarbutt Jul 16 '13 at 11:17
  • I've had problems doing intersections with flood layers before. The features are large and complicated. The way I've worked around it in the past is to split them into smaller features, so the spatial index can do more of the work. To do this, create a vector grid of the same extent as the flood layer (Vector > Research Tools > Vector Grid... Output grid as polygons), then intersect the grid with the flood layer. Then use the output instead of the flood layer in step 3. I'm guessing the reason the layer was empty is because it crashed. – Snorfalorpagus Jul 16 '13 at 18:00
  • Thanks again. The only problem now is that QGIS crashes every time I try to create the vector grid. I've followed the advice from here but it crashes every time. I've changed the parameters multiple times and tried only using the flood zone shapefile, instead of having my entire project file open and it fails every time. Any ideas? !Screenshot here. – KJGarbutt Jul 17 '13 at 11:26
  • The X and Y parameters you've specified are far too small. Try something like 1000 x 1000. You can even do this multiple times, i.e., do 5000 x 5000 first, use the output to create 500 x 500. See related answer here: gis.stackexchange.com/a/66319/12420 – Snorfalorpagus Jul 17 '13 at 15:24
  • I've nearly cracked it with your help! However, when I go to join the original LSOA layer with the intersected layer, I lose a lot of data. I think it is because some of the vector grid squares created fall within the same LSOA area and so have the same LSOA code as each each. Thus, I end up with 2+ percentage figures for each LSOA area when I do the join and I seem to only get one of them. Is there a way to sum each percentage for each vector grid square with the same LSOA? – KJGarbutt Jul 18 '13 at 14:02
6

You can use spatialite and some spatial SQL functions.

Select t1.geometry, t1.ID, area(t1.geometry), area(t2.geometry) ...... (anything you need to have in the table results)

(area(intersection(t1.geometry,t2.geometry))) as "Commun_AREA"

, ("Commun_AREA"*100/(area(t1.geometry))) as "Percent_AREA"

From lsoa as t1, flood_zone as t2

Where Intersects( t1.geometry,t2.geometry ) = 1
| improve this answer | |
3

This seems like something that could be done much easier than the answers submitted. I would use a simple python script personally:

floodName = "the layer name here"
boundryName = "the layer name here"
fieldName = "the name of the field to contain the output 1/0"
minCoverage = 0.5 # the minimum amount of area covered to write 1
updateMap = [] # this will store values to be written    

# get layers
floodLayer = QgsMapLayerRegistry.instance().mapLayersByName(floodName)[0]
boundryLayer = QgsMapLayerRegistry.instance().mapLayersByName(boundryName)[0]
fieldIndex = boundryLayer.dataProvider().fieldNameIndex(fieldName)    

# iterate through boundries
for b in boundryLayer.getFeatures():
    # get only flood features that intersect with this feature's bounding box
    # this will make the script go way faster than it would otherwise
    request = QgsFeatureRequest().setFilterRect(b.geometry().boundingBox())
    floodGeom = geometry()
    floodFeat = QgsFeature()
    iter = floodLayer.getFeatures(request)
    iter.nextFeature(feat)
    while iter.nextFeature(feat):
        floodGeom = floodGeom.combine(feat.geometry())
    intersectGeom = b.geometry().intersection(feat.geometry())
    if intersectGeom.area() > minCoverage * b.geometry().area():
        updateMap[b.id()] = {fieldIndex : 1}
    else:
        updateMap[b.id()] = {fieldIndex : 0}

boundryLayer.dataProvider().changeAttributeValues(updateMap)

this only evaluates flood polygons that intersect with the bounding box of each boundry layer so it should be fairly quick to run, then it only updates one field in the existing layer (instead of a complex operation of making a whole new layer and copying old values then deleting)

| improve this answer | |
2

I had the same problem as KJ following Snorfalorpagus' instructions using the "Intersect" method in Step 3. It took quite awhile to calculate and what I was left with was blank.

I tried following the same steps except using the "Clip" method in QGIS instead of Intersect -- so, in your example, what would be left would be the parts of the areas that are NOT covered by the flood zone. This seemed to work for some reason and I was able to use the "Area" field calculation from the previous step, plus a new "Area" calculation on the remaining parts of each polygon, to figure out the % of each area that was NOT covered by the other Polygon layer.

That's technically the reverse of what you asked for. But from there it's just a matter of subtracting each value from 1 to get what is covered by the flood zone.

| improve this answer | |
1

This question was posted some time ago, and the answer from @Snorfalorpagus was very thorough and effective (I used it several times!). I just wanted to note that as of QGIS 3.8, there is a new function called Overlap Analysis. This function:

"calculates the area and percentage cover by which features from an input layer are overlapped by features from a selection of overlay layers".

So, picking up off from step 2 in @Snorfalorpagus 's answer,

  1. Select Vector Analysis -> Overlap Analysis

  2. Input Layer (From example above): LSOA. Overlap Layer: Flood Areas.

This will create two new columns in LSOA:
- Column 1: actual area of overlap.
- Column 2: Percent area of overlap.

  1. You can easily convert to decimal (or proportion) from there using field calculator: (column_name / 100).
| improve this answer | |

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.