Issues with writing a dataframe to a shapefile

[ASandy]

Hello I've had an issue with writing a dataframe to a shapefile recently after being able to successfully complete the operation several times before. The error message I've now received is one I've never seen before.

This is the command that I have used before, now which is failing.

> #Write to a shapefile
> st_write(NYCHA_TotalCode, "NYCHA_TotalCode.shp")

This is the error message.
Error in CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  : 
  Layer creation failed.
In addition: Warning messages:
1: In abbreviate_shapefile_names(obj) :
  Field names abbreviated for ESRI Shapefile driver
2: In CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :
  GDAL Error 6: Geometry type of `3D Point' not supported in shapefiles.  Type can be overridden with a layer creation option of SHPT=POINT/ARC/POLYGON/MULTIPOINT/POINTZ/ARCZ/POLYGONZ/MULTIPOINTZ/MULTIPATCH.
> 
This csv file which I'm trying to turn into a shapefile has a geometry column among 6-7 others with 12,067 observations.
Please advise how I can remedy this as I've learned a lot in coding this year but still have a lot to learn regarding troubleshooting issues.

[technocrat]

Hi, it's hard to provide much specific help without a complete reproducible example, called a reprex. I do have a suggestion, however.

The sf (simple features) successor to the spatial package makes geometries much easier than dealing with shapefiles. In particular, they are data frames, with the geometry as one column. While sf can import properly formatted csv data, I don't see any direct export to shapefile function, so it may not work for you.

[ASandy]

OK, I understand. I will look into posting a decent reprex and into the sf package. Until then, can you hazard a guess as to why the "st_write" function worked fine for weeks before until yesterday and is now producing this error?

[jlacko]

I agree that we are missing a reprex, and there is not much help to be given without one.

Just an observation though: I have the impression @ASandy has been using sf::st_write().

This function does support writing in ESRI shapefile format, which can be advantageous when cooperating with backward colleagues who yet have to embrace R fully, and insist on working with legacy GIS platforms.

I cannot be certain this is the case, but I would not write the shp format off yet.

[ASandy]

Ok thank you for chiming in. I used the "st_read" function which I have been using on similar data for weeks. The "st_write" function used to work but I have included the relevant code to illustrate my issue. The error message is at the end. Please advise and let me know how I can you "help you help me"

getwd()

#Prepare all necessary or useful librairies
DaisyDuck <- c("sf","ggplot2","rstudioapi", "DescTools", "tibble",
               "readr", "ggspatial", "ggmap","rgdal", "rgeos","maptools",
               "GISTools", "dplyr", "tmap", "leaflet", "raster", "sp", "tidyverse")

#Load in all librairies at once
lapply(DaisyDuck, library, character.only = TRUE)

#Read in shapefile for complete study area by directory path
RF12067 <- st_read("H:\\VentilationUnits\\CombinedRoofFans\\ComboRF.shp")

#Begin cleaning data some of which is necessary from flawed GIS shapefile join
RF12067$TDS_NUM <- StrLeft(RF12067$LOCATION_C,3)

#Transforms coordinates of simple features dataframe based on certain projection and assigned to new object
RF12067_1 <- st_transform(RF12067, "+proj=longlat +ellps=WGS84 +datum=WGS84")

#Rearrange column order
FinalRF12067 <- RF12067_1[c("BOROUGH","DEVELOPMEN", "BLDG_NUM","TDS_NUM","LOCATION_C", "geometry")]

#Conversion of data type from factor to character for easier grouping
FinalRF12067$BLDG_NUM <- as.character(FinalRF12067$BLDG_NUM)
FinalRF12067$LOCATION_C <- as.character(FinalRF12067$LOCATION_C)
FinalRF12067$TDS_NUM <- as.character(FinalRF12067$TDS_NUM)

#Create a new column for each row representing an individual roof fan; it appends the tenant data system 
#code to a period and to the last 2 digits of the "Bldg_Num" 
FinalRF12067$NYCHA_BldgCode <-paste(FinalRF12067$TDS_NUM, sep = ".", StrRight(FinalRF12067$BLDG_NUM, 2))

#Create a table in which roof fan features are grouped by its' corresponding building code
LocationCount <- count(FinalRF12067, vars = NYCHA_BldgCode)

#Read column names and check datatype to prepare for cleaning and improved readability
names(LocationCount)


class(LocationCount$vars)


#Change the column names
names(LocationCount)[1] <- "NYCHA_BldgCode" 
names(LocationCount)[2] <- "No_RoofFans"

#Verify the datatype of each column
class(LocationCount$NYCHA_BldgCode)

class(LocationCount$No_RoofFans)

#For each column, concatenate with a comma between each element to prepare new dataframe in which 
#the roof fan counts are iterated individually to it's respective NYCHA property building
paste0(LocationCount$NYCHA_BldgCode, collapse = ",")
paste0(LocationCount$No_RoofFans, collapse = ",")

#Establish dataframe object with essential vectors from 'LocationCount' columns
LocationStringSet = data.frame(LocationCount$NYCHA_BldgCode, LocationCount$No_RoofFans)  

#Assign name to dataframe for function that takes columns of NYCHA building 
#and respective roof counts and creates a new list-columns
NYCHA_RF_LocationString = map2_df(LocationCount$NYCHA_BldgCode, LocationCount$No_RoofFans,
                                  ~data.frame(LocationCount.NYCHA_BldgCode=.x, LocationCount.No_RoofFans=.y, new.var=paste(.x, 1:.y, sep = ".RF")))

There were 50 or more warnings (use warnings() to see the first 50)

#Read out column names to begin improving readability
names(NYCHA_RF_LocationString)
names(NYCHA_RF_LocationString)[3] <- "RoofFanUniqueID"
names(NYCHA_RF_LocationString)
names(NYCHA_RF_LocationString)[1] <- "NYCHA_BldgCode"
names(NYCHA_RF_LocationString)[2] <- "RoofFanCt"

#Complete join that returns all rows of 1 table containing shapefile data 
#and the other table of corresponding NYCHA rooftop codes while filtering out enormous amount of duplicates
NYCHA_TotalCode <- merge(x= NYCHA_RF_LocationString, y= FinalRF12067[!duplicated (FinalRF12067$NYCHA_BldgCode),], by = "NYCHA_BldgCode", all.x= TRUE)

#Write to a shapefile
st_write(NYCHA_TotalCode, "NYCHArf_TotalCode.shp")

>Writing layer `NYCHArf_TotalCode' to data source `NYCHArf_TotalCode.shp' using driver `ESRI Shapefile'
Creating layer NYCHArf_TotalCode failed.
Error in CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  : 
                         Layer creation failed.
                       In addition: Warning messages:
                         1: In abbreviate_shapefile_names(obj) :
                         Field names abbreviated for ESRI Shapefile driver
                       2: In CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :
                                             GDAL Error 6: Geometry type of `3D Point' not supported in shapefiles.  Type can be overridden with a layer creation option of SHPT=POINT/ARC/POLYGON/MULTIPOINT/POINTZ/ARCZ/POLYGONZ/MULTIPOINTZ/MULTIPATCH.

[jlacko]

I suspect that there is an issue with your data. It somehow seems to have acquired a Z dimension (height). This is not supported by the shp format (it is positively ancient - based on dBASE file format, a leading database in the MS DOS era).

When testing the shapefile writing feature with a data frame of fake data (North Carolina is the Iris of spatial data) I get no funny messages, and the output is saved in point format (as expected).

library(sf)
library(tidyverse)

points <- data.frame(name = c("Raleigh", "Greensboro", "Wilmington"),
                     x = c(-78.633333, -79.819444, -77.912222),
                     y = c(35.766667, 36.08, 34.223333)) %>% 
  st_as_sf(coords = c("x", "y"), crs = 4326) %>% 
  st_transform(2264) # a feety CRS, just for the heck of it...

st_write(points, "sample_shapefile.shp")

[MartinStjernman]

I don't know if this is the problem but I note that you get a warning (number 2) which also give a GDAL Error 6: Geometry type of 3D Point not supported in shapefiles. I get a very similar error when trying to write to shapefile using an sf object that contains polygons with dimension XYZ (mysfobject$geometry[[1]][[1]][[1]] is a n by 3 matrix).
st_write(mysfobject, 'mysfobject.shp', delete_layer = TRUE)

Writing layer mysfobject' to data source mysfobject.shp' using driver ESRI Shapefile' Creating or updating layer test failed. Error in CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options), : Write error. In addition: Warning message: In CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options), : GDAL Error 6: Geometry type of 3D Multi Polygon' not supported in shapefiles. Type can be overridden with a layer creation option of SHPT=POINT/ARC/POLYGON/MULTIPOINT/POINTZ/ARCZ/POLYGONZ/MULTIPOINTZ/MULTIPATCH.

As my sfobject was created by reading in a shapefile obvioulsy having 3 dimensions using st_read() there is apparently not a problem to read such shapefiles, just writing them.
To fix it I have to remove the 3rd dimension by:
newsf <- st_zm(mysfobject, drop=T, what='ZM')
and then write to shapefile:
st_write(newsf, 'newsf.shp')

Writing layer newsf' to data source newsf.shp' using driver `ESRI Shapefile'
Writing 3978 features with 18 fields and geometry type Multi Polygon.

Hope this might be of some help.

/Martin