3

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?

watlbahamas plot result

Improve this question
8
  • 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
    Commented Jul 3, 2019 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
    Commented Jul 3, 2019 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
    Commented Jul 3, 2019 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
    Commented Jul 4, 2019 at 6:43
  • Which package does readSF come from?
    – Spacedman
    Commented Jul 4, 2019 at 6:44

1 Answer 1

Reset to default
8

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
Improve this answer
2
  • 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
    Commented Jul 5, 2019 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
    Commented Jul 5, 2019 at 22:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.