1

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?

  • If it were written, with all zero parameters it would have no effect at all. Would you like to test a principle or do you have some practical problem? – user30184 Aug 8 '18 at 18:51
  • See this @user30184 gis.stackexchange.com/questions/291734/… – Vijay Ramesh Aug 8 '18 at 19:18
  • There's a bunch of datums that could match +towgd84=0,0,0,0,0,0,0 and GRS80 and the US extent of the parameters: NAD83, NAD83 HARN, 2011, NSRS2007. If you don't look at the parameters, SIRGAS, NAD83 CSRS, etc. – mkennedy Aug 8 '18 at 20:07
  • @user30184 even with a +towgs84 parameter that isnt all zeroes I don't get it written to the shapefile .prj file. I'm wondering if its not part of the shapefile prj format... – Spacedman Aug 9 '18 at 6:43
  • @mkennedy wrote a few years ago in gis.stackexchange.com/questions/129764/… We don't support TOWGS84 nor some of the newer keywords. And that GDAL does not write it into .prj can be verified with command gdalsrsinfo -o wkt_esri "+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=1,1,1,1,1,1,1 +units=m +no_defs" – user30184 Aug 9 '18 at 8:16
1

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

  • Thanks @aldo_tapia, above worked i.e. adding +datum=NAD83. However, wondering if this will work for rasters and when using writeRaster(). I have to project a couple of rasters to this USGS albers equal area conic projection too and want to know if it will work. Do you know? – tg110 Aug 13 '18 at 0:07
  • @tg110 probably. This issue is writing vectors with OGR, rasters are written with GDAL – aldo_tapia Aug 13 '18 at 10:59
0

For some reason rgdal doesn't write the toWGS84 parameter to shapefiles:

If I create two identical objects in R, but with slightly different +towgs84 parameters:

> proj4string(d1)
[1] "+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=1,0,1,0,1,0,0 +units=m +no_defs"
> proj4string(d0)
[1] "+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"

then write to Shapefile, the .prj files are identical and show no trace of a datum shift parameter. ogrinfo tells me:

GEOGCS["GCS_GRS 1980(IUGG, 1980)",
    DATUM["unknown",
        SPHEROID["GRS80",6378137,298.257222101]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]],

for both files.

However, if I write them to a Geopackage,

> writeOGR(d0,"./d0.gpkg","d0","GPKG",overwrite=TRUE)
> writeOGR(d1,"./d1.gpkg","d1","GPKG",overwrite=TRUE)

I get the parameter stored and ogrinfo shows it:

For d0:

GEOGCS["GRS 1980(IUGG, 1980)",
    DATUM["unknown",
        SPHEROID["GRS80",6378137,298.257222101],
        TOWGS84[0,0,0,0,0,0,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]],

and for d1:

GEOGCS["GRS 1980(IUGG, 1980)",
    DATUM["unknown",
        SPHEROID["GRS80",6378137,298.257222101],
        TOWGS84[1,0,1,0,1,0,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]],

I'm not sure why this happens - I first thought that +toWGS84 datum shifts weren't supported in ESRI PRJ file format, but I have some shapefiles with TOWGS84 in them that seem to work fine:

DATUM["North_American_Datum_1927",
    SPHEROID["Clarke_1866",6378206.4,294.9786982138982,
        AUTHORITY["EPSG","7008"]],
    TOWGS84[-9.0,151.0,185.0,0.0,0.0,0.0,0.0],
    AUTHORITY["EPSG","6267"]],

But in any event, I'd recommend you use Geopackages from now on.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.