Convert shapefile into lat long points and export it as CSV in R?

Question

I have a shapefile that I extracted form a compilation of species distributions taken from IUCN.

> foo
Simple feature collection with 1 feature and 27 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -60.62262 ymin: -27.12379 xmax: -58.33061 ymax: -25.42827
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  id_no            binomial presence origin seasonal compiler year

I loaded the data with this

AR <- st_read("Mammals.shp")

foo <- AR %>% filter(binomial == "Foo bar")

I need to export these polygons, as points (lat and long) a .csv file to use the GeoCAT tools.

I tried this way but it didn't work,

foo_sf = st_as_sf(foo, coords = c("x", "y")) st_write(foo_sf, "foo.csv")

I am newbie to some of these packages so I don't know any workaround. I think the problem is the CSV needs lat and long points instead. I don't know if using SpatialPointsDataFrame would work here.

Answer 1

The help for st_write shows how to get the geometry for points into a CSV either as X,Y or as WKT, by using a specific layer_options parameter:

I have an sf object:

> dd
Simple feature collection with 1 feature and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.05609173 ymin: 0.642037 xmax: 0.05609173 ymax: 0.642037
epsg (SRID):    NA
proj4string:    NA
  z                    geometry
1 A POINT (0.05609173 0.642037)

then:

> st_write(dd, "out.csv", layer_options = "GEOMETRY=AS_XY")
Writing layer `out' to data source `out.csv' using driver `CSV'
options:        GEOMETRY=AS_XY 
features:       1
fields:         1
geometry type:  Point
> 

This produces the following CSV for my one-row data:

X,Y,z,
0.056091732578352,0.642037011450157,A

Now this is a for a POINT feature object - you can save a POLYGON feature object using WKT format, but is that what you want?

Normally to put a POLYGON of species occurrence into a model that needs points you'd take the centroid of the polygon, but if you really want to take the vertices of the polygon, then cast it into points before writing it as a CSV of points (as above):

Get one feature from the NC dataset:

> nc = st_read(system.file("shape/nc.shp", package="sf"))
> nc = nc[1,]

This is a polygon with some attributes:

> nc
Simple feature collection with 1 feature and 14 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
epsg (SRID):    4267
proj4string:    +proj=longlat +datum=NAD27 +no_defs
   AREA PERIMETER CNTY_ CNTY_ID NAME  FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
1 0.114     1.442  1825    1825 Ashe 37009  37009        5  1091     1      10
  BIR79 SID79 NWBIR79                       geometry
1  1364     0      19 MULTIPOLYGON (((-81.47276 3...

Casting to POINT results in a number of features with repeated attributes:

> st_cast(nc, "POINT")
Simple feature collection with 27 features and 14 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
epsg (SRID):    4267
proj4string:    +proj=longlat +datum=NAD27 +no_defs
First 10 features:
    AREA PERIMETER CNTY_ CNTY_ID NAME  FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
1  0.114     1.442  1825    1825 Ashe 37009  37009        5  1091     1      10
2  0.114     1.442  1825    1825 Ashe 37009  37009        5  1091     1      10
3  0.114     1.442  1825    1825 Ashe 37009  37009        5  1091     1      10
4  0.114     1.442  1825    1825 Ashe 37009  37009        5  1091     1      10
5  0.114     1.442  1825    1825 Ashe 37009  37009        5  1091     1      10
 [etc]

Then:

> pts = st_cast(nc,"POINT")
Warning message:
In st_cast.sf(nc, "POINT") :
  repeating attributes for all sub-geometries for which they may not be constant
> st_write(pts,"/tmp/test_pts.csv",layer_options="GEOMETRY=AS_XY")

will get you:

X,Y,AREA,PERIMETER,CNTY_,CNTY_ID,NAME,FIPS,FIPSNO,CRESS_ID,BIR74,SID74,NWBIR74,BIR79,SID79,NWBIR79
    -81.4727554321289,36.2343559265137,0.114,1.442,1825,1825,Ashe,37009,37009,5,1091,1,10,1364,0,19
    -81.5408401489258,36.2725067138672,0.114,1.442,1825,1825,Ashe,37009,37009,5,1091,1,10,1364,0,19
    -81.5619812011719,36.2735939025879,0.114,1.442,1825,1825,Ashe,37009,37009,5,1091,1,10,1364,0,19
    -81.6330642700195,36.3406867980957,0.114,1.442,1825,1825,Ashe,37009,37009,5,1091,1,10,1364,0,19
    -81.7410736083984,36.3917846679688,0.114,1.442,1825,1825,Ashe,37009,37009,5,1091,1,10,1364,0,19

Now I think the tool you are using needs columns called "Latitude" and "Longitude" instead of X an Y (read the docs in detail for the file format). I'm not sure if this can be specified in st_write or the CSV OGR driver, so you may have to edit the file.

Answer 2

Try as.data.frame(st_coordinates(x)), which will return a data.frame object of [x,y] coordinates, and then write.csv() You can also pull a column from the sf object following something like data.frame(ID = x$ID, st_coordinates(x))

The sp SpatialPointsDataFrame class approach, using rgdal to read the shapefile, would follow:

library(sp)
library(rgdal)
foo <- readOGR(getwd(), "myshp")
write.csv( data.frame(ID=foo$ID, coordinates(foo)), "myshp.csv", row.names = FALSE )