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 a decently detailed shapefile with polygon/multipolygon features (the file is about 500mb). It's actually a shapefile of the entire world, with the features representing coastlines. I need to divide this data using a grid. To be clear, I don't want to 'sort' the data, but actually cut the polygons up into tiles. I realize this question has been asked before but the solutions I found didn't work for me.

I've tried:

  • Using QGIS and intersecting my shapefile contents with a vector grid -- the results are terrible. Most of the major landmass magically disappears, though it seems like smaller chunks of land sometimes make it. I should note that this method works really well with much simpler data (ie. less points)

  • Using OGR's Intersection tools. I tried it both through ogr2ogr and even by rolling my own C++ tool. Both of them have the same problem as QGIS. They also don't exhibit this problem for simple files, but fail the more complex ones. For reference, I'm using a shapefile of Australia and New Zealand, under 20mb in size, and both QGIS and OGR fail to 'gridify' it.

Someone suggested using PostGIS at one point, since it has an intersection function -- but PostGIS's ST_Intersect uses the same GEOS backend as OGR does. In fact they both call the same function as far as I can tell, so I don't think that PostGIS will yield different results.

I was looking for suggestions as to what else I could try. I need a robust application or toolkit that can divide highly detailed shapefiles into tiles.

EDIT: Adding some more information

In response to Simbamangu:

  • The shapefile is basically coastline data from OpenStreetMaps. It's a merged version of the 'processed_p' file (so its not split up into tiles) that I got by emailing their dev list. Note that their splitting up of tiles (into 100km x 100km chunks with overlap) isn't necessarily what I want -- I don't want overlap, and I want the freedom to choose the grid size, or I'd just use the default processed_p.

  • By default the coastline data has geometry errors reported by QGIS. I fix these errors with a small tool I put together using some code I found designed to specifically address this problem (repairing geometry errors in coastline data: https://github.com/tudelft-gist/prepair). Running over the files with this tool fixes virtually all the errors QGIS picks up. I only attempt to do the intersection after cleaning the files.

  • Exactly what I did using QGIS: Open the data to make sure it looks fine in QGIS. Try dividing it into tiles by creating a layer of tiles using Vector Grid with a specified spacing, and then intersecting the two layers -- no go. Try using a smaller data set -- select features in Oceania (Aus, NZ) to try a smaller data set -- this shape file is < 20mb in size. Again try dividing it, doesn't work.

  • What I did with OGR: ogr2ogr directly using the '-spat' and '-clipsrc' options with spat_extent. Also wrote a small C++ tool that works on WKT, so I convert the shapefile to WKT using ogr2ogr, then feed the text file to my application. It runs through the file and calls the Intersection() method documented here: http://www.gdal.org/ogr/classOGRGeometry.html. I think it ends up doing the exact same thing as using ogr2ogr directly.

In response to Brent:

  1. It does. Everything is in WGS84 Lat/Lon
  2. I would have thought that the opposite is true -- that for a given set of grid tiles, it would take way longer to intersect one giant multipolygon rather than a bunch of fragmented features which could be more spatially localized to each tile, but this is an interesting suggestion -- I'll try it and report back.
  3. No attribute fields are kept during the process, I'm only interested in geometry.
  4. I'm not sure, but I think you're saying I should select the polygons that overlap a given grid tile and then perform the intersection. This is too cumbersome manually with QGIS. My tool already does this to a certain extent with a bounding box check. There's a bit of a speed up, but the end result is still poor and not noticeably different.
  5. This isn't an option. Right now I'm trying to divide the data up so that its 1 deg lat x 1 deg lon, and I'm looking for a general/robust methodology that works with all cases. I've tried increasing the grid size (ie 10x10) to see if I'd get better results and I don't see any correlation between grid size and the quality of the output.

Edit #2:

I've tried playing around with this more and in general it just seems that the results are unreliable both using GEOS and with QGIS (which uses fTools, I don't know if that in turn uses GEOS again). I was wrong in stating the size of the grid has nothing to do with the results -- the bigger the grid is, the better the results (that's good to know but still not a solution). Here's a screenshot of a really spaced out grid that mostly worked, but failed partially in one tile:

enter image description here

The geometry is clean -- QGIS shows 0 errors with the "Check Validity" tool. I'm not looking to approach this problem in a step by step fashion; verifying whether or not certain features failed the intersection on a dataset this big when its not visually apparent (and it won't be with smaller tiles) isn't practical.

share|improve this question
    
Where did you get the world or Australia shapefile? I'd suspect that the geometry of that file may have some issues (try Vector|Geometry Tools|Check Geometry Validity in QGIS). Have just tried an intersect on a smaller world shapefile and 5 degree tiles and it works perfectly in QGIS. –  Simbamangu May 14 '12 at 4:42
1  
Tried this with the 100K Australian coastline from Geoscience Australia (20MB) and 4 degree tiles, also works fine (QGIS 1.7.4, OSX 10.7). Could you describe in more detail your data and what you did? –  Simbamangu May 14 '12 at 6:51
    
Thanks for all the extra information. I suspect there's something odd about the OSM data; try it with the dataset I mentioned and see if you get better results. I seem to remember experiencing some weirdness with OSM lake data in the past, will try to look it up. –  Simbamangu May 14 '12 at 16:36
    
Could you share the dataset, or even a clipped part of it (like in your example above)? –  Simbamangu May 17 '12 at 11:54
    
mediafire.com/?196nrdlabcek4xk –  Pris May 17 '12 at 18:58
add comment

3 Answers

It definitely sounds like you have geometry issues. It is unlikely that will be able yield clean results from a dirty input file regardless of the software used, unless you first address your geometry problems. Once you get your geometry issues sorted out, you could try the following if you are still having issues:

1) Make sure that your grid dataset has the same projection as your world polygon dataset. If not, recreate it in the proper projection.

2) Convert all features to single part - much easier to process

3) Remove all extraneous fields keeping only id field which will enable you to join your attributes back after the intersection has been performed - again much easier to process

4) Instead of intersecting the entire grid dataset with the entire world polygon dataset, try loop over your grid polygons, selecting the intersecting polygons in your world dataset and the performing a clip to based on your grid polygon. This will enable you to isolate any problems and in the end you can merge results together to achieve your original goal.

5) Try using larger grid polygons.

share|improve this answer
    
+1 Really interesting - how much does it affect geoprocessing speed if you keep the ID field, or multipart, in the data? –  Simbamangu May 14 '12 at 15:56
1  
I've never actually tried to quantify the differences. I can only speak from experience where overaly geoprocesing operations failed and these are the sort of things that helped solve the problem. –  Brent Edwards May 14 '12 at 16:45
    
I was unsuccessful in getting (2) to work at all. Selecting features and trying to merge them using QGIS basically seems to lock up my system -- maybe its still processing stuff, but at that rate its not practical: I left my system on overnight with QGIS still trying to merge a couple of features in the dataset and it was still going at it in the morning. –  Pris May 16 '12 at 18:27
1  
There should not be any merging involved. The goal is to explode multipart features. For example, in your screenshot of the failed tile, the goal is to explode all of your records containing grouped, spatially disjoint polygons like the island features along the coast of BC and Alaska, into separate, single part polygon records. This can be achieved in QGIS using the "Multipart to singleparts" tool under the Vector>Geometry Tools menu. –  Brent Edwards May 17 '12 at 12:53
    
Once your convert to single part feature, you should re-validate your geometry, just to be certain that everything is clean. –  Brent Edwards May 17 '12 at 14:01
show 1 more comment
up vote 4 down vote accepted

I just ended up creating my own tools to do this.

I used the Clipper library (http://www.angusj.com/delphi/clipper.php) along with OGR to divide my data set up. Something to note is performing intersections naively with this lib takes very long, so I instead used a quadtree approach... ie, divide into four grid cells, divide each of those into four more, etc, until you get your desired resolution. The lib works great though, I've attached a screenshot showing the results on the eastern hemisphere:

enter image description here

The above result took about 4.5 hours on a 1.33GHz processor.

Here are the tools in case someone runs into a similar issue in the future. Please note that they're hacked together proof-of-concepts and you probably shouldn't use them directly (might serve as a good starting point for something though): https://github.com/canurabus/scratch/tree/master/polytk

share|improve this answer
    
The linked code is no longer available :-( –  Shaun McDonald Apr 16 at 13:10
add comment

Another approach might have been to try a vector-to-raster conversion to create a point dataset & then use the point dataset as a basis for writing some code to create your tiles.

Seems like you've solved the problem though. Just my two cents worth

Ian

share|improve this answer
add comment

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.