- I have been using RGDAL to transform and write out spatial data in GDA94
- It is a problem to know what I need to do to create the right prj file for ArcGIS to read without complaining
- I have an example of code below, that I used on a dataset I knew was in GDA94 when I read it in. I want to do this mostly for when I have to convert from one to another, but I have done the manual hack a couple of times now and thought I better just check with an expert.
R Code:
infile <- "ap_map"
outfile <- gsub("map", "map_GDA94.shp", infile)
outfile
setwd(indir)
shp <- readOGR(".", infile)
setwd(projdir)
plot(shp, add = T)
#str(shp)
shp@proj4string@projargs
#[1] "+proj=longlat +ellps=GRS80 +no_defs"
# confirm this is GDA94
epsg <- make_EPSG()
str(epsg)
epsg[grep(4283, epsg$code),]
# code note prj4
# 212 4283 # GDA94 +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
# Arcmap sees this as GRS 1980(IUGG, 1980) which I think is the same thing
# to be on the safe side I will force it to refer to GDA94
shp2 <- spTransform(shp, CRS("+init=epsg:4283"))
shp2@proj4string@projargs
#now write
setwd(outdir)
dir()
writeOGR(shp2, outfile, gsub(".shp", "", outfile), driver = "ESRI Shapefile")
setwd(projdir)
# Checking this shows it did not write the correct prj file, but I believe that this is because the GDA94 definition is not different to the WGS80 params.
# one other option is to manually replace that prj file with the correct text found at http://spatialreference.org/ref/epsg/4283/
# SO I did this to avoid any future confusions
# OLD GEOGCS["GRS 1980(IUGG, 1980)",DATUM["D_unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
# NEW GEOGCS["GDA94",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
# checked with ArcMap and they align ok
I asked on of my colleagues and here is his reply
- I usually use the sp and maptools library, rather than the readOGR and writeOGR functions. Generally my workflow, as an example of projecting a GDA94 zone 55 file to WGS84, would be something like:
R Code:
shp <- readShapePoly("myfile.shp")
proj4string(shp) <- "+proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
shp2 <- spTransform(shp, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))
writePolyShape(shp2,"myfile_p.shp",proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))
-
I’ve never noticed any problems doing this, the projection file is generated and I’ve never noticed any alignment problems in ArcMap. It’s true that I recall ArcMap sometimes displays an incorrect plaintext projection description, but things seem to be in the right place.
-
Also, I don’t really ever rely on EPSG numbers in R - I just grab the proj4 strings and try to use them directly.
TODO, see if PostGIS handles this OK
- I wonder if my FOSS GIS stack should still pass things through PostGIS as I suspect it does these things better.
- I might check the output if I run it through the DB use st_transform and then extract to shapefile.
