Creating Shapefiles in R

Question

I'm trying to create a shapefile in R that I will later import to either Fusion Table or some other GIS application.

To start,I imported a blank shapefile containing all the census tracts in Canada. I have attached other data (in tabular format) to the shapefile based on the unique ID of the CTs, and I have mapped my results. At the moment, I only need the ones in Vancouver and I would like to export a shapefile that contains only the Vancouver CTs as well as my newly attached attribute data.

Here is my code (some parts omitted due to privacy reasons):

shape <- readShapePoly('C:/TEST/blank_ct.shp') #Load blank shapefile
shape@data = data.frame(shape@data, data2[match(shape@data$CTUID, data2$CTUID),]) #data2 is my created attributes that I'm attaching to blank file
shape1 <-shape[shape$CMAUID == 933,] #selecting the Vancouver CTs

I've seen other examples using this: writePolyShape to create the shapefile. I tried it, and it worked to an extent. It created the .shp, .dbf, and .shx files. I'm missing the .prj file and I'm not sure how to go about creating it. Are there better methods out there for creating shapefiles?

Any help on this matter would be greatly appreciated.

This is how I write my files `writeOGR(obj = opno.skupaj.mean[[1]], dsn = "q:/path/to/file/spat_odstrel_skupaj.shp", layer = "spat_odstrel_skupaj", driver = "ESRI Shapefile")`. Notice that the layer and filename are identical (minus the `.shp`). — Roman Luštrik, Apr 12 '13 at 07:45

mnel · Answer 1 · 2013-04-12T07:48:43.013

8

Use rgdal and writeOGR. rgdal will preserve the projection information

something like

library(rdgal)

shape <- readOGR(dsn = 'C:/TEST', layer = 'blank_ct')
# do your processing
shape@data = data.frame(shape@data, data2[match(shape@data$CTUID, data2$CTUID),]) #data2 is my      created attributes that I'm attaching to blank file
shape1 <-shape[shape$CMAUID == 933,]
writeOGR(shape1, dsn = 'C:/TEST', layer ='newstuff', driver = 'ESRI Shapefile')

Note that the dsn is the folder containing the .shp file, and the layer is the name of the shapefile without the .shp extension. It will read (readOGR) and write (writeOGR) all the component files (.dbf, .shp, .prj etc)

edited Apr 12 '13 at 07:48

answered Apr 12 '13 at 03:16

mnel

113,303
27
265
254

1

That won't work as-is, you've not specified the `shape` or `shape1` object in the `writeOGR` call - how does it know what to write? – Spacedman Apr 12 '13 at 07:16
1

I think the dsn path to `writeOGR` needs to be a full file path ending in `.shp`. – Roman Luštrik Apr 12 '13 at 07:46
@Roman it does not need that, can be writeOGR(x, ".", "layername", "ESRI Shapefile") or writeOGR(x, "C:/temp/layername.shp", "layername", "ESRI Shapefile") – mdsumner Apr 12 '13 at 11:35
@mdsumner ha! I wonder what I did wrong. I was unable to write with your method, but now it works. Did we cross the zone where strange things happen... very often? – Roman Luštrik Apr 12 '13 at 16:41
1

@mdsumner It's important that the path in `dsn` does not include a trailing slash if one decides to not specify the full file path. – Roman Luštrik Apr 13 '13 at 16:39

score 4 · Answer 2 · answered Apr 12 '13 at 14:54

Problem solved! Thank you again for those who help!

Here is what I ended up doing:

As Mnel wrote, this line will create the shapefile.

writeOGR(shape1, dsn = 'C:/TEST', layer ='newstuff', driver = 'ESRI Shapefile')

However, when I ran this line, it came back with this error:

Can't convert columns of class: AsIs; column names: ct2,mprop,mlot,mliv

This is because my attribute data was not numeric, but were characters. Luckily, my attribute data is all numbers so I ran transform() to fix this problem.

shape2 <-shape1
shape2@data <- transform(shape1@data, ct2 = as.numeric(ct2),
mprop = as.numeric(mprop),
mlot = as.numeric(mlot),
mliv = as.numeric(mliv))

I tried the writeOGR() command again, but I still didn't get the .prj file that I was looking for. The problem was I didn't specified the coordinate systems for the shapefile when I was importing the file. Since I already know what the coordinate system is, all I had to do was define it when importing.

readShapePoly('C:/TEST/blank_ct.shp',proj4string=CRS("+proj=longlat +datum=WGS84")

After that, I re-ran all the things I wanted to do with the shapefile, and the writeOGR line for exporting. And that's it!

Creating Shapefiles in R

2 Answers2

Linked