2

I'm trying to use the workflow described by @aflaciano at How to render area that crosses 180°? I did the following:

  1. downloaded and installed gdal-111-1600-x64-core.msi from here on my Windows 7 64 bit pc.
  2. downloaded data from http://www.naturalearthdata.com
  3. used the following batch-file adapted from @aflaciano:
    rem Clip world.shp between -169° and +180° of longitude in order to exclude the Chukotka Peninsula (1st part):
    rem ==========================================================================================================
    ogr2ogr world_clip.shp ne_110m_admin_0_countries.shp -clipsrc -169 -90 180 90
    rem Shift and clip world.shp between +180° and +191° of longitude in order to include the Chukotka Peninsula (2nd part):
    rem ====================================================================================================================
    ogr2ogr Chukotka_Peninsula.shp ne_110m_admin_0_countries.shp -dialect sqlite -sql "SELECT sov_a3, ShiftCoords(geometry,360,0) FROM ne_110m_admin_0_countries" -clipsrc 180 -90 191 90
    rem Merge the two parts:
    rem ====================
    ogr2ogr world_shifted.shp world_clip.shp
    ogr2ogr -update -append world_shifted.shp Chukotka_Peninsula.shp -nln world_shifted
    rem Finally, dissolve countries boundaries of world_shifted.shp obtaining world_output.shp:
    rem =======================================================================================
    ogr2ogr world_output.shp world_shifted.shp -dialect sqlite -sql "select ST_Union(Geometry),sov_a3 FROM world_shifted GROUP BY sov_a3"

At step 3 last line of the batch-file I must have done an error, because dissolve does not work for me, but I get no error message:

ogr2ogr world_output.shp world_shifted.shp -dialect sqlite -sql "select ST_Union(Geometry),sov_a3 FROM world_shifted GROUP BY sov_a3"

I tried also other variables for dissolve, e.g. sovereignt or admin, geounit, adm0_a3

enter image description here

Glorfindel
  • 1,096
  • 2
  • 9
  • 14
Kurt
  • 7,087
  • 4
  • 33
  • 51

1 Answers1

3

This was a nice problem. Fortunately you gave excellent step-by-step instructions how to repeat the issue.

The last step with ST_Union does not fail and therefore you do not get any error message. The union is created but the function does not "dissolve" the parts but writes them as separate parts of a multipolygon. And that happens because the common line segment only looks like common. Closer look at the coordinates reveals that they are not equal. This is the left side:

180 64.97970870219838,  
180 68.96363636363661,

And this is the right side:

180.00000000000003 64.97970870219837, 
180.00000000000003 68.96363636363637,

Because of the small rounding errors ST_Union writes the left side and right side into separate members of a multipolygon. I do not know why it fails for both of us but not for @aflaciano. Perhaps we have different platforms or GDAL and Spatialite versions.

I made some further trials with good and bad news. Good news is that it is possible to combine the halves into one polygon with Spatialite by using a ST_SnapToGrid function. Bad news is that at least my GDAL dowloaded from gisinternals.com is not compiled with enough recent Spatialite and SnapToGrid is missing from it.

What should be done is to round the coordinates before executing ST_Union so that vertices match exactly. The command to use would be then

ogr2ogr world_output.shp world_shifted.shp -dialect sqlite -sql "select ST_Union(ST_SnapToGrid(Geometry,0.00001)),CNTRY_NAME FROM world_shifted GROUP BY CNTRY_NAME"

GDAL trunk version from gisinternals downloaded today (June 2nd, 2014) throws an error

no such function: ST_SnapToGrid

However, the command works with spatialite-gui v. 1.7.1 and the result is good.

Chukotka Peninsula though Spatialite and ST_Union + ST_SnapToGrid

user30184
  • 65,331
  • 4
  • 65
  • 118
  • thank you for your kindly help and explaining, why my approch fails. Could you please give more details howto do the job with saptialite-gui v. 17.1. I did the following:
    1. created new db and imported shapefile world_shifted.shp
    2. run your query select ST_Union(ST_SnapToGrid(Geometry,0.00001)),CNTRY_NAME FROM world_shifted GROUP BY CNTRY_NAME

    (I used sov-a3 instead of CNTRY_NAME, btw) 3) no errors, less rows than before query, so it seems ok question of a rockie: how can I save this query into a new db or shapefile, so that I can use it in software like qgis ???

     – Kurt Jun 02 '14 at 16:49
  • This worked for me. 1) create an empty Spatialite db 2) import shapefile "world_shifted" 3) create a new table with SQL create table output as select CastToMultiPolygon(ST_Union(ST_SnapToGrid(Geometry,0.00001))) as geometry, sov_a3 FROM world_shifted GROUP BY sov_a3 4) refresh the table of contents in spatialite-gui, select column "geometry" from the output table with right click, run Recover geometry column (4326, XY, multipolygon) 5) Refresh. Now you can export the table as shapefile with right click on "geometry" or open directly with QGIS. – user30184 Jun 02 '14 at 17:15
  • following your steps: step 1: create empty spatialite db- ok step 2: import world_shifted.shp -ok step 3: doing your query create table output as select CastToMultiPolygon(ST_Union(ST_SnapToGrid(Geometry,0.00001))) as geometry, sov_a3 FROM world_shifted GROUP BY sov_a3 gives the following error:

    SQL query returned an empty ResultSet This is not an error

    any hints?

     – Kurt Jun 02 '14 at 17:34
  • It's not an error. Query did not return any resultset because it created a new table. Refresh with right click on top of TOC and you will see the new table on the list. – user30184 Jun 02 '14 at 17:43
  • now after step 4 (recover geometry column (4326,XY,multipolygon) I get the following error: Geometry column doesn't satisfies required constraints a ROLLBACK was automatically performed :-( – Kurt Jun 02 '14 at 17:54
  • doing UPDATE "output" SET "geometry" = SetSrid("geometry", 4326) WHERE Srid("geometry") = 0 allowed me to do your step 4 afterwards successfully. But: after exporting to shapefile and importing into qgis I dont have the parts as separate parts of a multipolygon, no I have them as different polygons, that means just like in the file world-shifted :-( . maybe in SnapToGrid(Geometry,0.00001)) I should increase the value?? @user30184 which file from naturalearth did you use (I use ne_110m_admin_0_countries), so that I can retry again – Kurt Jun 02 '14 at 18:11
  • correction to my last post: after doing UPDATE "output" SET "geometry" = SetSrid("geometry", 4326) WHERE Srid("geometry") = 0 and doing step 5 (= refresh) I could export the shapefile, which is now correctly dissolved :-) – Kurt Jun 02 '14 at 18:42
  • Great. If you had set the SRID when you imported the shapefile you would have avoided this SetSrid step. And with manual geometrytype selection "MULTIPOLYGON" the CastToMultipolygon step. – user30184 Jun 02 '14 at 18:48