I am trying to create a shapefile from a .csv file in R. I can successfully create a working shapefile with this code, but when I try to view the attribute table for it in arcMap 10.4, the table only contains columns for objectID, latitude, and longitude. I need the info in the other columns for symbolizing the data correctly.
Is there a way to specify which columns of the data to keep in the shapefile?
(I know how to make a data frame into a shapefile within arcMap, but I would like to automate this step in R if possible.)
library(maptools)
library(rgdal)
library(sp)
#data
site <- c("a","b","c","d")
prop_c <- c(0.88,0.48,0.15,0.47)
prop_b <- c(0.17,0.18,0.09,0.08)
minus_c <- 1-prop_c
minus_b <- 1-prop_b
lat <- c(44.22,38.38,33.35,43.48)
long <- c(-124.45, -123.70, -124.40, -124.05)
MyData <- cbind.data.frame(site, prop_c, prop_b, minus_c, minus_b, lat, long)
#convert data to shapefile
WGScoor<- MyData #data to convert
coordinates(WGScoor)=~long+lat #column names of the lat long cols
proj4string(WGScoor)<- CRS("++proj=longlat +datum=WGS84") # set coordinate system to WGS
WGScoor.df <- SpatialPointsDataFrame(WGScoor, data.frame(id=1:length(WGScoor)))
LLcoor<-spTransform(WGScoor.df,CRS("+proj=longlat"))
LLcoor.df=SpatialPointsDataFrame(LLcoor, data.frame(id=1:length(LLcoor)))
writeOGR(LLcoor.df, dsn=getwd(),layer="MyShapefile",driver="ESRI Shapefile")
#successfully creates shapefile that opens in arcGIS,
#but shapefile does not include other columns in original data frame
2 Answers 2
You are overcomplicating things.
> WGScoor<- MyData
> coordinates(WGScoor)=~long+lat
> proj4string(WGScoor)<- CRS("+proj=longlat +datum=WGS84")
at this point you have something you can save as a shapefile with all the columns intact, but you seem to want to project out of WGS84, so let's do that:
> LLcoor<-spTransform(WGScoor,CRS("+proj=longlat"))
and let's save this:
> raster::shapefile(LLcoor, "MyShapefile.shp")
I'm using here the shapefile
function from the raster
package. Its just a nice wrapper to the writeOGR
function you used, but makes things a bit easier.
Now let's see what columns are in the resulting shapefile via the command line:
$ ogrinfo -so -al MyShapefile.shp
INFO: Open of `MyShapefile.shp'
using driver `ESRI Shapefile' successful.
Layer name: MyShapefile
Geometry: Point
Feature Count: 4
Extent: (-124.450000, 33.350000) - (-123.700000, 44.220000)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
DATUM["unknown",
SPHEROID["WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]]
site: String (80.0)
prop_c: Real (24.15)
prop_b: Real (24.15)
minus_c: Real (24.15)
minus_b: Real (24.15)
which shows they are all there. I'm not sure your transform from "+proj=longlat +datum=WGS84"
(and you had an extra +
in your string) to "+proj=longlat"
actually does anything, but I've left it.
-
Thank you! That really helps clear up the steps of the code, and now the shapefile saves with the columns intact.Mon Mo– Mon Mo2016年10月14日 23:10:36 +00:00Commented Oct 14, 2016 at 23:10
-
@Spacedman What is the explanation of this
CRS("+proj=longlat +datum=WGS84")
? I cannot understandMK Huda– MK Huda2022年11月07日 00:50:53 +00:00Commented Nov 7, 2022 at 0:50 -
It uses the
sp
package functions to convert a string"+proj=longlat +datum=WGS84"
to a coordinate reference system object.CRS("+proj=longlat +datum=WGS84")
- note this Q is 6 years old and if you are coding new things now try and use thesf
package. If you still don't understand something try posting as a new question.Spacedman– Spacedman2022年11月07日 08:18:35 +00:00Commented Nov 7, 2022 at 8:18
When you use the function SpatialPointsDataFrame(), the "data.frame" argument is where you get all the data you want in your shape. Now, you're putting only an ID.
Try:
LLcoor.df2=SpatialPointsDataFrame(LLcoor, data.frame(id=1:length(LLcoor), MyData))
You can also do later on:
LLcoor.df2@data <- MyData
However this is risky as you can get order problems. Make sure the other of your shape and your data frame is teh same before replacing it that way.
Hope it helps,
-
This method works too. It creates a few additional columns in the attribute table when opened in ArcGIS, but all the essential columns are there.Mon Mo– Mon Mo2016年10月14日 23:14:31 +00:00Commented Oct 14, 2016 at 23:14