Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have two polygon shapefiles, each of which contain multiple (possibly overlapping) polygons.

Is there some way, using GRASS, to calculate the number of polygons belonging to shapefile A that fall within each of the polygons that exist in shapefile B? Preferable would be a method that doesn't require rasterization, but I'm open to going down that route if necessary.

(I'm aware that this can be done in R, using the over() function, but these are large files I'm referring to, and so I'd prefer to do it in GRASS.)

share|improve this question
add comment

2 Answers

up vote 4 down vote accepted

To continue webrian's suggestion:

v.select ain=poly_a b_in=poly_b oper=within out=a_within_b

Now the output vector "a_within_b" will contain those polygons from a that were totally within b. But be careful: when you import a shapefile into GRASS with overlapping polygons, the overlaps will be split up by the topology "clean" operation, so you'll have more, smaller polygons in the GRASS vector than in the original shapefile.

Now to count the number of polygons, you might (depending on the database GRASS is using) be able to do:

echo 'SELECT count(cat) FROM a_within_b' | db.execute

But this raises the question: why not do this right in a spatial database from the start? Import both shapefiles into spatialite, and run:

SELECT count(poly_a.pk_uid) FROM poly_a, poly_b 
WHERE ST_Within(poly_a.Geometry, poly_b.Geometry);
share|improve this answer
    
I tried out spatialite and it seems like a viable option. I was interested in identifying the features of A (or their count) that occurred within (incl. crossing, etc.) each of the features of B, I ended up using: SELECT polyA.pk_uid AS pk_A, polyB.pk_uid AS pk_B FROM polyA, polyB WHERE ST_Intersects(polyB.Geometry, polyA.Geometry) AND polyB.ROWID IN ( SELECT pkid FROM idx_polyB_Geometry WHERE pkid MATCH RTreeIntersects ( MBRminX(polyA.geometry), MBRminY(polyA.geometry), MBRmaxX(polyA.geometry), MBRmaxY(polyA.geometry) )); Comments? –  jbaums Jun 28 '12 at 23:42
    
The above still seems quite slow (although maybe my expectations are too high...). Can anyone see a way to speed this up? –  jbaums Jun 28 '12 at 23:51
    
If you're using a recent version of spatialite, then the format for adding a spatial index to the "WHERE" clause is simpler: ...AND polyB.ROWID IN (SELECT ROWID FROM SpatialIndex WHERE f_table_name='polyB' AND search_frame=polyA.Geometry); Maybe that will speed things up? –  Micha Jun 29 '12 at 21:40
add comment

Have a look at the v.select module in GRASS GIS. There is no need to rasterize.

share|improve this answer
    
Thanks webrian, but I can't quite figure out how to use it for my purposes... I get the feeling it just returns those features in ainput that overlap with features of binput (overlap is the default, and apparently the only option available without GEOS, which I haven't yet worked out how to use with GRASS via R). I will try directly from GRASS though, to see whether there might be more to it than I realise. Can you provide an additional info re how to invoke v.select? –  jbaums Jun 26 '12 at 12:03
add comment

We're looking for long answers that provide some explanation and context. Don't just give a one-line answer; explain why your answer is right, ideally with citations. Answers that don't include explanations may be removed.

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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