[R] Re ading and Creating Shape Files

cls59 chuck at sharpsteen.net
Thu Oct 29 15:28:02 CET 2009



PDXRugger wrote:
> 
> Hello R Community,
>    I have imported a dataset which contain X Y coordinates and would like
> to recreate a shape file after some data analysis.  What i have done is to
> import some taxlot data and join them based on some criteria.  I want to
> check to see how well the joining went by reviewing the results in GIS.  
> 
> A couple things.  I cant seem to import a shape file correctly using the
> maptools package and the readShapeSpatial.  I have tried 
> 
> Building=file("data/input/BuildingShape/Building.shp")
> Bldg<-readShapeSpatial(fn=data/input/BuildingShape/Building,proj4string=NAD83)
> #----------------------
> Bldg<-readShapeSpatial(data/input/BuildingShape/Building,proj4string=NAD83)
> #---------------
> Building=file("data/input/BuildingShape/Building.shp")
> Bldg<-readShapeSpatial(Building,proj4string=NAD83)
> 
> I know i am mis interpreting the documentation but it doesnt seem like it
> is very complicated so i am of course confused.
> 
> 

I haven't used the maptools package for this kind of operation, so I offer
any specific advice. However, the value of proj4string seems to be a little
odd-- unless NAD83 is a variable that contains a string. It should probably
be a quoted list of PROJ4 declarations such as:

proj4string = '+proj=longlat +datum=NAD83'

If maptools is leveraging the sp package, then you probably need to enclose
the string in the CRS() function:

proj4string = CRS('+proj=longlat +datum=NAD83')

Note that the above example assumes coordinates are in lat/lon using the
NAD83 datum. If your data is in a different projection, such as UTM, you
will need to change the PROJ4 string accordingly.



PDXRugger wrote:
> 
> 
> Also, i am wondering if i can create a shape file by simply using XY
> coordinates from a data frame.
> So for:
> 
> 
> Ycoord=c( 865296.4, 865151.5, 865457.0 ,865363.4 ,865311.0, 865260.9
> ,865210.7 ,865173.3,
> 865123.6 ,865038.2 ,864841.1 ,864745.4 ,864429.1 ,864795.6 ,864334.9
> ,864882.0)
> 
> Xcoord=c( 4227640 ,4227816 ,4228929 ,4228508 ,4229569 ,4229498 ,4226747,
> 4226781, 4229597,
> 4229204, 4228910, 4228959 ,4229465 ,4229794 ,4229596 ,4229082)
> 
> Lot<-c(1900 , 2000,  2100  , 100   ,200  , 300,   400 ,  500 ,  600 ,  701
> ,  900 , 1000 , 1100,
>   300   ,100,   200)   
> 
> XYcoord<-spCbind(Ycoord,Xcoord) #doesnt work so
> 
> XYcoord=c(Ycoord,Xcoord)
> 
> TaxLots<-cbind(Ycoord,Xcoord,Lot)
> 
> writeSpatialShape(XYcoord, TaxLots..,
> file=data/input/test/Taxlots,strictFilename=FALSE)
> 
> 
> 
> So help reading in shape files and then creating them using XY coordinates
> if possible
> Any help would be appreciated.  Thank you.
> 
> 
> 

Maybe maptools provides a nice way to do this-- again I haven't used it
much. My shapefile workflow usually centers on the sp and rgdal packages.
First, sp is used to create a spatial object that holds the coordinates and
data. Assuming you have point data:

  require(rgdal)
  lots <- SpatialPointsDataFrame( coords = cbind(Xcoord,Ycoord), data =
data.frame( Lot = Lot ))

*Note that SpatialPointsDataFrame also takes a proj4string argument of the
form:

proj4string = CRS( 'proj declarations' )

I have omitted it since I don't know what projection your data is in.


You can then create a shapefile using the writeOGR() routine in rgdal:

  writeOGR( lots, dsn = 'tstShapefile', layer = 'tstShapefile', driver='ESRI
Shapefile')


The readOGR() function can also be used to read a shapefile-- note that you
give it the name of the directory containing the shapefile components and
not the name of an individual component such as 'shapefile.shp'.

Hope this helps!

-Charlie

-----
Charlie Sharpsteen
Undergraduate
Environmental Resources Engineering
Humboldt State University
-- 
View this message in context: http://www.nabble.com/Reading-and-Creating-Shape-Files-tp26098828p26114143.html
Sent from the R help mailing list archive at Nabble.com.




More information about the R-help mailing list