2

Trying to polygonize large rasters* but since the process is taking forever/crashing in GDAL, I am resorting to doing so in ArcGIS. However there is a 2GB limit on shapefiles in Arc so I am taking a combination approach where I split the rasters into 4 parts with gdalwarp, polygonize them in ArcGIS (polygonizing the quarters in GDAL is also taking forever while the process was extremely fast in Arc), and then putting them back together in GDAL to get around the slowness of gdal_polygonize.py and also the 2GB ArcGIS limit.

I am basing my solution partially on the answer found here, but am running into a problem with the final Dissolve step.

*The pixels of the rasters I am polygonzing only contain 1's (i.e. the pixel value of the things I'm interested in as polygons) and NoData (do not wish to have polygons in these areas), so the attribute table after polygonization looks like this:

attribute table

Which is want I expect.

The next step is to merge these (4) polygons back together in GDAL like this:

ogr2ogr -f "ESRI Shapefile" merge.shp seg7_9__h06v01_1_2011to2015_min_a.arcpoly.shp -nln merge # create merge.shp
ogr2ogr -f "ESRI Shapefile" -update -append merge.shp seg7_9__h06v01_1_2011to2015_min_b.arcpoly.shp -nln merge # add b.shp to merge.shp
ogr2ogr -f "ESRI Shapefile" -update -append merge.shp seg7_9__h06v01_1_2011to2015_min_c.arcpoly.shp -nln merge
ogr2ogr -f "ESRI Shapefile" -update -append merge.shp seg7_9__h06v01_1_2011to2015_min_d.arcpoly.shp -nln merge

Which works fine and gives me a result that looks like this, where the highlighted polygon is from one of the input shapefile parts and shows the edge between the two:

merge shp

My last step/goal is to dissolve the merge.shp so that the polygons on the edge that are touching (such as the one highlighted in red and the portion above it) become one polygon. In the post linked above, they suggested using the following command (command1):

ogr2ogr "output.shp" "input.shp" -dialect sqlite -sql "SELECT ST_Union(geometry), field FROM input GROUP BY field"

Because I do not really understand SQL or SQLITE, I wanted to try this command (command2) at first just to see what would happen:

ogr2ogr out_dissolve1.shp merge.shp -dialect sqlite -sql "SELECT ST_Union(geometry), field FROM merge"

but got the following error:

ERROR 1: In ExecuteSQL(): sqlite3_prepare(SELECT ST_Union(geometry), field FROM merge):
 wrong number of arguments to function ST_Union()

I looked online and tried looking at documentation of ST_Union but could not find anything useful.

In conclusion:

  1. Can anyone explain what the SQL statement of the function from command1 referenced above does and

  2. Does anyone know how I can adjust the ogr2ogr command to obtain the desired result, which is dissolving the touching/intersecting polygons into one?

UPDATE: A couple of other variations I've tried and the errors:

I thought using the GRIDCODE field would do the same thing I want because every polygon has a GRIDCODE=1

ogr2ogr dissolve.shp merge.shp -dialect sqlite -sql "SELECT ST_Union(geometry), GRIDCODE FROM merge GROUP BY GRIDCODE"

ERROR 1: In ExecuteSQL(): sqlite3_prepare(SELECT ST_Union(geometry), GRIDCODE FROM merge GROUP BY GRIDCODE):
  wrong number of arguments to function ST_Union()

and

ogr2ogr dissolve.shp merge.shp -dialect sqlite -sql "SELECT ST_Union(geometry) FROM merge" 

ERROR 1: In ExecuteSQL(): sqlite3_prepare(SELECT ST_Union(geometry) FROM merge):
  wrong number of arguments to function ST_Union()

Not sure why I'm getting this error or how to fix it.

UPDATE2: I think I found the right command to dissolve based solely on geometry, but I'm getting an error with the ST_Union command:

ogr2ogr dissolve.shp merge.shp -dialect sqlite3 -sql "SELECT ST_Union(geometry) as geometry FROM merge"

ERROR 1: Undefined function 'ST_Union' used.

I tried using GUnion and got the same error.

  • instead of shapefiles, look into feature classes, no size limit on those. – Maksim May 5 '17 at 19:47
  • You talk about rasterizing but are you rather polygonizing? – user30184 May 5 '17 at 19:49
  • @user30184 cr*p you are right, I am trying to polygonize but keep getting mixed up. thanks. edited. – user20408 May 5 '17 at 19:51
  • How much is large as rows and columns? – user30184 May 5 '17 at 20:42
  • @user30184 The original rasters are 55500x55500 before splitting into quarters. The polygonization in arc of these split rasters is veryyyy fast though – user20408 May 5 '17 at 20:46
3

If I understand correctly, the problem crops up when you try to Dissolve the results shapefile which you get after combining 4 smaller shapefiles via either Append or Merge.

Possibility 1: Update. Rather than use Append or Merge on those 4 files, use Update from the Analysis Tools\Overlay Toolbox. Set the keep_borders variable to NO_BORDERS and those lines should disappear.

Possibility 2: How are you at python scripting? If you've any familiarity with any programming language, you'll be able to do decent at it, even with no training. Also, ModelBuilder includes the option to export a model as a python script. The reason I ask is that a python script doesn't have ESRI's unchangeable 2GB limit. If you have a step or steps that should work, but don't, a pure python script using ArcPy should provide a workaround.

Possibility 3: You might be able to work around this error by running the dissolve command after each shapefile is Appended to the results file. I honestly don't know if that'd work, but ArcGIS is often best approached in a "This shouldn't work, but I'll try it anyway, just in case" state of mind.

Possibility 4: Take the 4 files, create layers from each, use select by location to find only those shapes that overlap or intersect each other, then dissolve just those features. That might take enough of a load off your memory for the process to complete. You'll now have a results shapefile. Use the same select by location (just invert the selection) to select all the features that didn't get dissolved, and Append them to that same results file.

Possibility 5: Eliminate. Merges smaller polygons into the larger ones that they border/intersect.

Possibility 6: Repair Geometry. Those borders might not trip the self-intersection detector in Check Geometry. But, if they do, Repair Geometry might fix it.

I hope one of these ideas works.

  • Hi, Thanks so much for your thoughtful answer. However I tried several of these methods and no matter what, ArcGIS would not create a a SHP greater than 2GB, so it looks like I'm definitely going to have to do this in GDAL/OGR, unfortunately. – user20408 May 5 '17 at 23:18
  • I meant to add in my comment that even running this in python outside of arc (with arcpy) yields these same results. It does not give me an error during processing but stops making the shp past 2GB and when you add the result to ArcMap a "could not draw" error does arise. – user20408 May 5 '17 at 23:32
  • Have you tried converting the SHP to a feature class in a file GDB? – P.T. Curran May 8 '17 at 21:35
  • I have not tried this, no. What's the best way to approach this if the merged shapefiles I'm trying to dissolve are larger than 2GB? Can I merge the individual inputs into a GDB then dissolve? Thank you – user20408 May 9 '17 at 14:26
  • Nevermind. I tried merging the inputs into a GDB and it worked, however when I tried to Dissolve I got an "Out of memory" error. Any idea on how to get around that? Thanks. – user20408 May 9 '17 at 16:16

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.