I am writing out a shp, in ESRI shapefile driver. In R, the projection of this shp is
+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0
+y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
I wrote this shp out as below-
writeOGR(x_projected, dsn="address path","name of shp to be stored", driver="ESRI Shapefile", overwrite_layer=TRUE)
However, on reading this projected, written shp back in, I see that the datum specified by +towgs84=0,0,0,0,0,0,0 does not exist. Also, when plotting this shp in ArcMap or QGIS, the datum is known confirming that the datum has not been written out.
Why is writeOGR() not writing out the datum i.e. why is the +towgd84=0,0,0,0,0,0,0 component missing?
Answer
Check the following wiki:
Datum
You can use the +nadgrids-Parameter to select an existing gridfile or define your own datum with the +towgs84 Parameter to use a 3- or 7-Parameter Transformation. To use the buildin "datums" (or dati?) you can use the +datum-Parameter and select one of the DatumCode in the following list from pj_datums.c (e.g. +datum=NAD27).
Is weird, defining +towgs84=
in R, as:
library(sp)
library(rgdal)
data(meuse)
coordinates(meuse) <- ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")
meuse_nad83 <- spTransform(meuse, "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ")
writeOGR(meuse_nad83, dsn="/path/to","name", driver="ESRI Shapefile", overwrite_layer=TRUE)
You got this from QGIS:
But, if you define +datum=
:
meuse_nad83b <- spTransform(meuse, "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ")
writeOGR(meuse_nad83b, dsn="/to/path","name", driver="ESRI Shapefile", overwrite_layer=TRUE)
You got:
I can't answer why writeOGR()
doesn't work, just recommend you to use +datum=
parameter
No comments:
Post a Comment