R help - Nov 2011 - How to write a shapefile with projection

Hi,

?

I am trying to write a shapefile with projection. I have
my data in a data.frame called try and consists in xy coordinates and a
numerical attribute value z1.

?

Libraries loaded are: sp, rgdal, raster, maptools

?

head(try)

???????? x?????? ??? y??? ? ? ? ? ? ? ??? z1

1 610237.1 ???????? 3375751
???????????? 8.221

2 610236.1 ???????? 3375750
???????????? 8.153

3 610236.1 ???????? 3375749
???????????? 8.275

4 610236.1 ???????? 3375748
???????????? 8.251

5 610236.1 ???????? 3375747
???????????? 8.217

6 610236.1 ???????? 3375746
???????????? 8.196

?

#Get the projection from a raster ?named llev I have loaded before:

crs <- projection(llev)

?

# get a spatial point data frame from my data

crest.sp <- SpatialPointsDataFrame(try[,1:2], try,
proj4string=CRS(crs))

?

summary(crest.sp)

Object of class SpatialPointsDataFrame

Coordinates:

??????? min?????? ??????? max

x? 610235.1? ?????? 610354.1

y 3374862.4 ?????? 3375751.4

Is projected: TRUE 

proj4string :

[+proj=utm +zone=15 +ellps=GRS80 +datum=NAD83 +units=m
+no_defs +towgs84=0,0,0]

Number of points: 890

Data attributes:

????? ???????? ?x????????? y???????????????? z1?????? 

?Min.?? :? 610235?? Min. : 3374862?? Min.: 6.966? 

?1st
Qu.:610269?? 1st Qu.:3375085?? 1st Qu.:7.570? 

?Median
:610298?? Median :3375307?? Median :7.901? 

?Mean?? :610300??
Mean? ?:3375307??
Mean?? :7.882? 

?3rd
Qu.:610334?? 3rd Qu.:3375529?? 3rd Qu.:8.180? 

?Max.?? :610354??
Max.?? :3375751?? Max.??
:8.756? 

?

#write the shapefile

writePointsShape(crest.sp, "I:/LA_levee/Shape/llev_crest_pts6")

?

If I load this shapefile in ArcGIS is has no projection.
I also looked at the write.shapefile command from shapefiles library and again I
get a file without projection. Is there any way to write the projection for the
shapefile in R? Probably I can do something in ArcGIS, but I would like to have
my shapefile from R complete.

?

Thanks,

?

Monica

?

> Hi,
>
> Sorry i have put such a detailed question to the list about writing a
shapefile with projection. I realized that if i use writeOGR from rgdal and not
the other write shapefile functions i can get a shapefile with projection
recognized by ArcGIS. The command is (in case anybody wonders):
>
> ?writeOGR(crest.sp, "I:\\LA_levee\\Shape",
"llev_crest_pts6", driver = "ESRI Shapefile")
>
> where crest.sp is a spatial point data frame with projection.
>
> Thanks,
>
> Monica
Indeed.

writePointsShape() does not write the projection file, but  using the 
function showWKT from rgdal, you can also create one like that:

writePointsShape(crest.sp,"crest")
cat(showWKT(proj4string(crest.sp)),file="crest.prj")

Patrick

Hi Patrick,

Thanks for letting me know. I mostly use rgdal to read and write rasters so
until now i kind of ignore other functionality. Unfortunately i supposed that a
package dedicated to shapefiles would be the answer and had the functionality i
needed. But rgdal does a nice job in saving my files as i need. It is good to
know how to add the projection file to the shapefiles for the future, if it is
not generated from the onset.

Thanks again,

Monica

----------------------------------------
> Date: Sat, 5 Nov 2011 14:27:29 +0100
> From: patrick.giraudoux at univ-fcomte.fr
> To: pisicandru at hotmail.com
> CC: r-help at r-project.org
> Subject: re: How to write a shapefile with projection
>
> > Hi,
> >
> > Sorry i have put such a detailed question to the list about writing a
shapefile with projection. I realized that if i use writeOGR from rgdal and not
the other write shapefile functions i can get a shapefile with projection
recognized by ArcGIS. The command is (in case anybody wonders):
> >
> > ?writeOGR(crest.sp, "I:\\LA_levee\\Shape",
"llev_crest_pts6", driver = "ESRI Shapefile")
> >
> > where crest.sp is a spatial point data frame with projection.
> >
> > Thanks,
> >
> > Monica
>
> Indeed.
>
> writePointsShape() does not write the projection file, but using the
> function showWKT from rgdal, you can also create one like that:
>
> writePointsShape(crest.sp,"crest")
> cat(showWKT(proj4string(crest.sp)),file="crest.prj")
>
> Patrick
>
>