0

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.

  • As far as I can tell from the FAQ kew.org/geocat-faqs#2 the tool can only import point data, so why are you trying to import polygons from a shapefile? – Spacedman Sep 4 '18 at 17:47
  • Yes, @Spacedman I need to convert the polygon to points first. I have the species distribution as a polygon and I need to convert it. I need points lat and long as you noted, do you know how to do this? – Andrés Parada Sep 4 '18 at 17:52
  • But the tool on the web site only takes species locations as points. If you have a polygon for the species and convert it to points that's like saying your animals are only found on the corners of the polygon, Is that right? – Spacedman Sep 4 '18 at 18:48
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.

  • Nice shortcut with st_write and layer_options = "GEOMETRY=AS_XY" – Jeffrey Evans Sep 4 '18 at 16:31
  • This code exports a csv with shape_Leng and shape_Area 6.2788635563 2.44299406288 and not points as I need, additionally I got this message > Warning message: In CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options), : GDAL Message 1: Geometry type Multi Polygon is not compatible with GEOMETRY=AS_XYZ. – Andrés Parada Sep 4 '18 at 16:58
0

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 )

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.