1

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 http://www.gisinternals.com/sdk/ 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

|improve this question
2

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

|improve this answer
  • 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 2 '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 2 '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 2 '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 2 '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 2 '14 at 17:54

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.