2

I'm trying to rasterize the "Decreasing land productivity" shapefile from https://wad.jrc.ec.europa.eu/geoportal, using:

gdal_rasterize -3d -co COMPRESS=DEFLATE -ts 43200 21600 ./lpd_int2.shp ./ldp_int2.tif

(also tried without the COMPRESS=DEFLATE)
It has been running for 6 days and is only 40% completed.

Though the output is 43200x21600, this still seems like an excessively slow process. Is there some obvious thing I'm doing which is making it take too long?

$ ogrinfo -ro -al -listmdd lpd_int2.shp
INFO: Open of `lpd_int2.shp'
      using driver `ESRI Shapefile' successful.

Layer name: lpd_int2
Metadata domains:
  (default)
Metadata:
  DBF_DATE_LAST_UPDATE=1919-09-29
Geometry: Polygon
Feature Count: 226
Extent: (-179.999989, -55.600000) - (179.991667, 72.400000)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433],
    AXIS["Longitude",EAST],
    AXIS["Latitude",NORTH]]
Data axis to CRS axis mapping: 1,2
id_country: Real (33.15)
name_iso31: String (254.0)
iso3: String (254.0)
sqkm: Real (33.15)
lpd_sqkm: Real (33.15)
lpd_p: Real (33.15)
shape_leng: Real (33.15)
shape_area: Real (33.15)
OGRFeature(lpd_int2):0
  id_country (Real) = 3.000000000000000
  name_iso31 (String) = Afghanistan
  iso3 (String) = AFG
  sqkm (Real) = 643490.775693999952637
  lpd_sqkm (Real) = 52103.569999999999709
  lpd_p (Real) = 8.100000000000000
  shape_leng (Real) = 90294575.479599997401237
  shape_area (Real) = 80173734283.800003051757812
  MULTIPOLYGON  ...etc, etc...
  • The shapefile seems to contain extremely heavy geometries - only 226 polygons but the size of the .shp part of shapefile is 467 MB. – user30184 Oct 25 '19 at 20:14
  • One of the monster multipolygons contains 516716 rings! – user30184 Oct 25 '19 at 20:30
  • Oh. That seems like a lot. I'll start searching for how to simplify shapefile geometry. – DGentry Oct 25 '19 at 21:06
2

The source dataset has huge multipolygons and GDAL has difficulties to deal with them. That applies also to ogr2ogr that cannot explode multipolygons into single part polygons with option -explodecollections. However, I found a reasonable simple workaround.

  1. Convert the shapefile into SpatiaLite database with ogr2ogr

    ogr2ogr -f sqlite -dsco spatialite=yes lpd_int2.sqlite lpd_int2.shp -nlt promote_to_multi

  2. Explode the multipolygons into single part polygons with spatialite-gui. SpatiaLite has a special SQL function for the purpose as documented in http://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html

    Will create a new out_table directly corresponding to in_table. The output table will be arranged in such a way that each row will always contain an elementary Geometry; so each time that a MULTI-type Geometry is found in the input table it will be split into mamy distinct rows. out_pk is the name to be set for the output Primary Key, and out_multi_id is the name to be set for a second column within the output table where to store the original Primary Key. The optional argument transaction determines if an internal SQL Transaction should be automatically started or not (the default setting if not explicitly overridden is TRUE).

    Syntax to use:

    SELECT ElementaryGeometries('lpd_int2','GEOMETRY','exploded','id_1','id_2');

    It took 2 min 30 seconds to create the new table with my laptop.

  3. Run gdal_rasterize by using the exploded polygons from SpatiaLite as source data.

    gdal_rasterize -burn 255 -co tiled=yes -ts 43200 21600 lpd_int2.sqlite -l exploded lpd_int2.tif

    I modified your command a bit because the source shapefile nor the table of exploded polygons is not 3D and there are no Z coordinates to use as burn values. It took about 15 minutes to create a GeoTIFF and it looks like this in QGIS

enter image description here

EDIT:

The issue with ogr2ogr and -explodecollections is now also fixed in GDAL master, see https://github.com/OSGeo/gdal/issues/1959.

  • 1
    Thank you very much for this. I don't have a Windows system to run spatialite-gui, but was able to produce the same result using: ogrinfo ./lpd_int2.sqlite -sql "SELECT ElementaryGeometries('lpd_int2','GEOMETRY','exploded','id_1','id_2');" – DGentry Nov 2 '19 at 14:10

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.