Saturday 21 December 2019

coordinate system - writeOGR() not writing out datum of shp


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:


enter image description here


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:


enter image description here


I can't answer why writeOGR() doesn't work, just recommend you to use +datum= parameter


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...