0

I'm building map viewer using d3js and data from WFS service.
In Poland we have 3 levels of administrative divisions:

  • province
  • district
  • community

Using WFS service I can get exact boundaries of each level. For example for provinces I'm using http://sdi.geoportal.gov.pl/wfs_prg/wfservice.aspx?request=getFeature&version=1.1.0&service=WFS&typename=gmgml:WOJEWODZTWA
For districts: http://sdi.geoportal.gov.pl/wfs_prg/wfservice.aspx?request=getFeature&version=1.1.0&service=WFS&typename=gmgml:POWIATY

Then using ogr2ogr I'm converting those GML files to TopoJSON and displaying using d3js. This part works fine.

Right now I need to know which district belongs to which province, but unfortunately this isn't available in file that comes from WFS.
I need this to be able to filter my map by provinces - display only districts for specific provinces.

I know that I can open SHP that file in QGIS and add attribute for every district, but I would have to add this to 2000+ items (before that checking manually if that district belongs to that province), and I would like to avoid that extra work.

Can this be done automatically? Using some kind of software?
I have only tiny experience with QGIS so I didn't found that option.

i imagine is as checking if calculated center point o district belongs to province boundaries, but maybe this is simpler task?

2

What you want to do is create a spatial join in QGIS.

Load your two shape files into QGIS. Then select Vector->Data Management Tools->Join Attributes by Location. Select your district layer as the Target Vector Layer, your provinces layer as Join Vector Layer. For Attribute Summary select "Take Attributes of first located feature". Select your output layer. For Output Table, select "Only Keep Matching Records". If the polygons are accurate, it shouldn't matter which value you select for Attribute Summary or Output Table.

The result of this operation will be a new layer with as many rows as your district layer. Each row will contain all the fields of the district layer PLUS the fields of the Province layer. You can then delete any unneeded fields that were inherited from the Province table (presumably everything except the province name or id).

Your first sanity check is that the output layer should have exactly as many rows as your input (Districts) layer. The more difficult check is that the district is in the right province. If the boundaries of the two layers don't quite align, a district might be seen as being in two different provinces, and the join will select the first one seen. One way to address this to add a new integer field to the provinces layer. Call it something like "Verify". Set this field to 1 for every record. When doing the join, set the Attribute Summary to "Take Summary of Intersecting Features" and select "Sum". Then when you're checking the resulting join, look for records where the Verify field is greater than 1. These are the records you need to check and perhaps fix manually.

An alternative approach that directly deals with districts intersecting multiple provinces is by selecting the province with the greatest area in the intersection. The best I could come up with is to use Vector->GeoProcessing->Intersect to intersect all the districts with all the provinces, then take the ratio of the area of the intersection to the area of the district. (This requires you have an area field in each of the two input tables.) There's a sharp cutoff between those that are nearly 1 and those that are nearly 0. Delete the ones that are nearly zero, which should leave a table with districts and provinces. Note this table is NOT the correct answer yet - since these districts are the result of the intersection with provinces, so it's not the same as your district input file. I would then delete all fields in the intersection table except the district ID and the province ID, then do an ordinary join on the district ID field with your original district table.

(Sorry for this somewhat rambling answer, but I'm trying to figure out how to do this as I go.)

  • Thanks for such detailed answer, unfortunately data that comes from this government service isn't accurate enough and I can't get any good results. I didn't check second method, but I'll give it a try tomorrow. My idea for now is to check if centroid of district belongs to province. First I must check if centroid of each district is located inside district borders and then do a loop to check if that point belongs to province shape. This probably isn't build in, so I'll have to learn python to script this. Thanks again for grate answer, it helped me to know QGIS more :) – Misiu Jan 12 '14 at 20:13

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.