st_write fails for object created with st_difference using R?

Question

Using publicly available sea area shapefiles

bahamas: http://www.marineregions.org/gazetteer.php?p=details&id=8404 ("eez.zip")
watl: http://www.marineregions.org/gazetteer.php?p=details&id=7135 ("iho_quadrants.zip")
gom: http://www.marineregions.org/gazetteer.php?p=details&id=4288

I joined the Bahamas to the Western Atlantic

watl <- read_sf("iho_quadrants_20150810.shp")
bahamas <- read_sf("eez.shp")
watlbah <- st_union(watl, bahamas)

Then removed overlapping areas with the Gulf of Mexico

watlbahamas <- st_difference(watlbah, GOM)

and the Slope Sea. File: https://drive.google.com/open?id=1eDmn8k3cHqhrqdLdLr8wKDYmnDRid_NQ

slopesea <- read_sf(slopesea.GPKG)
watlbahamas <- st_difference(watlbahamas, slopesea)

Then saved the result

st_write(obj = watlbahamas, dsn = "watlbahamas.GPKG")

Which resulted in the following:

Error in st_write.sf(obj = watlbahamas, dsn = "watlbahamas.GPKG") : 
  failed writing to temporary file /tmp/RtmpgG9ykj/file15ca4a4711fe.GPKG
In addition: Warning messages:
1: In CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :
  GDAL Error 1: sqlite3_exec(CREATE TABLE "watlbahamas" ( "fid" INTEGER PRIMARY KEY AUTOINCREMENT NOT NULL, "geom" GEOMETRY, "name" TEXT, "id" TEXT, "bearings" TEXT, "name_1" TEXT, "mrgid_1" MEDIUMINT, "name_2" TEXT, "mrgid_2" MEDIUMINT, "name_3" TEXT, "mrgid_3" MEDIUMINT, "mrgid" REAL, "geoname" TEXT, "mrgid_ter1" REAL, "pol_type" TEXT, "mrgid_sov1" REAL, "territory1" TEXT, "iso_ter1" TEXT, "sovereign1" TEXT, "mrgid_ter2" REAL, "mrgid_sov2" REAL, "territory2" TEXT, "iso_ter2" TEXT, "sovereign2" TEXT, "mrgid_ter3" REAL, "mrgid_sov3" REAL, "territory3" TEXT, "iso_ter3" TEXT, "sovereign3" TEXT, "x_1" REAL, "y_1" REAL, "mrgid_eez" REAL, "area_km2" REAL, "name.1" TEXT, "id.1" TEXT, "longitude" REAL, "latitude" REAL, "min_x" REAL, "min_y" REAL, "max_x" REAL, "max_y" REAL, "area" REAL, "mrgid.1" MEDIUMINT, "Id" MEDIUMINT)) failed: duplicate column name: Id
2: In CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :
  GDAL Error 1: failed to prepare SQL: SELECT "fid", ST_MinX("geom"), ST_MaxX("geom"), ST_MinY("geom"), ST_MaxY("geom") FROM "watlbahamas" WHERE "geom" NOT NULL AND NOT ST_IsEmpty("geom")
3: In CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :
  GDAL Error 1: sqlite3_exec(DROP TABLE "watlbahamas") failed: no such table: watlbahamas
4: In CPL_write_ogr(obj, tmp, layer, driver, as.character(dataset_options),  :
  GDAL Error 1: sqlite3_exec(CREATE TABLE "watlbahamas" ( "fid" INTEGER PRIMARY KEY AUTOINCREMENT NOT NULL, "geom" GEOMETRY, "name" TEXT, "id" TEXT, "bearings" TEXT, "name_1" TEXT, "mrgid_1" MEDIUMINT, "name_2" TEXT, "mrgid_2" MEDIUMINT, "name_3" TEXT, "mrgid_3" MEDIUMINT, "mrgid" REAL, "geoname" TEXT, "mrgid_ter1" REAL, "pol_type" TEXT, "mrgid_sov1" REAL, "territory1" TEXT, "iso_ter1" TEXT, "sovereign1" TEXT, "mrgid_ter2" REAL, "mrgid_sov2" REAL, "territory2" TEXT, "iso_ter2" TEXT, "sovereign2" TEXT, "mrgid_ter3" REAL, "mrgid_sov3" REAL, "territory3" TEXT, "iso_ter3" TEXT, "sovereign3" TEXT, "x_1" REAL, "y_1" REAL, "mrgid_eez" REAL, "area_km2" REAL, "name.1" TEXT, "id.1" TEXT, "longitude" REAL, "latitude" REAL, "min_x" REAL, "min_y" REAL, "max_x" REAL, "max_y" REAL, "area" REAL, "mrgid.1" MEDIUMINT, "Id" MEDIUMINT)) failed: duplicate column name: Id
5: In CPL_write_ogr(obj, tmp, layer, driver, as.character(dataset_options),  :
  GDAL Error 1: failed to prepare SQL: SELECT "fid", ST_MinX("geom"), ST_MaxX("geom"), ST_MinY("geom"), ST_MaxY("geom") FROM "watlbahamas" WHERE "geom" NOT NULL AND NOT ST_IsEmpty("geom")
6: In CPL_write_ogr(obj, tmp, layer, driver, as.character(dataset_options),  :
  GDAL Error 1: sqlite3_exec(DROP TABLE "watlbahamas") failed: no such table: watlbahamas

So I tried to simplify the object (https://gis.stackexchange.com/a/243576/16071):

watlbahamas <- rmapshaper::ms_simplify(input = as(watlbahamas, 'Spatial')) %>%
  st_as_sf()

But the install fails (I can debug this but want to know if I'm on the right track first). Trying a different way:

watlbahamas <- rgeos::gSimplify(watlbahamas, tol = 0.1) %>%
  st_as_sf()
# Error in rgeos::gSimplify(watlbahamas, tol = 0.1): no slot of name "proj4string" for this object of class "sf"

I feel like I've done a relatively simple operation but either there are elements which cause downstream problems in st_difference, or, more likely, I've gone about this the wrong way.

So (notwithstanding some typos from my trying to remove some text bulk) it seems the slope sea file is the culprit, and that st_difference on only GOM works fine.

Could it be that two+ st_difference calls create >1 set of the same named columns?

You've got lots of missing " marks. You've written st_write(obj = watlbahamasTMP, dsn = watlbahamas.GPKG) but not shown us what watlbahamasTMP is, and you've not quoted the file name. Saving watlbahamas as in your code works for me so I think you've probably not really done what you say you've done.... — Spacedman, Jul 03 '19 at 19:21
Sorry, "TMP" was left in in error, was just a quick temporary rename attempt while things were failing. Filename was quoted originally but as part of a paste0() with folder path that I figured would just bog down the text unnecessarily, again now corrected. Waiting on results of another code run then will retry. — dez93_2000, Jul 03 '19 at 19:51
Ah rats, so it seems the slope sea poly is what's killing it - I figured since it was the same operation it wouldn't make a difference. Question now edited. — dez93_2000, Jul 03 '19 at 20:10
Also you don't define GOM but I assume its GOM = st_read("iho.shp") where iho.shp is the download from the third of your list of links? — Spacedman, Jul 04 '19 at 06:43
re GOM, correct, sorry for the omission. re read_sf, it's from sf, an alias of st_read; looks like I should probably use that instead. — dez93_2000, Jul 05 '19 at 21:47
ah good spot, that was again me omitting lines for brevity then manually adding them back in later; corrected — dez93_2000, Jul 07 '19 at 16:41

score 7 · Accepted Answer · answered Jul 04 '19 at 06:58

7

GPKG column names must be unique and are case insensitive. Do you have any columns where the case-folded names are the same?

> any(duplicated(toupper(names(watlbahamas))))
[1] TRUE

Yes. Which ones?

> which(duplicated(toupper(names(watlbahamas))))
[1] 42
> names(watlbahamas)[42]
[1] "Id"
>

So name 42 is a duplicate of something, and a quick inspection shows that's a duplicate (apart from upper/lower case) of name number two, "id". Let's change that one:

> names(watlbahamas)[2] 
[1] "id"

> names(watlbahamas)[2] = "id2"

And test....

> st_write(watlbahamas,"wid.gpkg")
Writing layer `wid' to data source `wid.gpkg' using driver `GPKG'
features:       2
fields:         42
geometry type:  Multi Polygon

And it worked.

If you read the error message you'll see there's a massive clue:

[...stuff...] failed: duplicate column name: Id

answered Jul 04 '19 at 06:58

Spacedman

63,755
5
81
115

Good stuff, thanks. So as I suspected, doing st_difference >1X will create dupe id colnames that make st_write fail. I was hoping a simplify operation would flatten various features and thus address the dupe colname issue but a simple rename is wiser. I wonder if this is something that's worth raising with Edzer & the sf team... – dez93_2000 Jul 05 '19 at 21:49
No, your spatial data doesn't have duplicated column names because R is case-sensitive for column named - you have "Id" and "id" and R treats them differently. But some spatial data formats may not allow this. At the R level, the code doesn't know if the output format driver allows it or not, and doesn't care. Its job is to pass the data to the output format driver and if the driver raises any error messages then its the users job to read those to figure out what went wrong. The R level code can't know the individual restrictions of all the output drivers. – Spacedman Jul 05 '19 at 22:09

st_write fails for object created with st_difference using R?

1 Answers1