[R-sig-Geo] +no_uoff flag in proj4string, epsg 26931
Roger Bivand
Roger.Bivand at nhh.no
Tue Aug 12 21:40:08 CEST 2014
On Tue, 12 Aug 2014, Gregovich, Dave P (DFG) wrote:
> Hello,
> I would like to project data in epsg 26931 (NAD83 Alaska State Plane
> zone 1 meters), variant A. There are two variants to the projection:
> --variant A has its origin on the central line of the projection, on the
> 'aposphere equator'
> --variant B has its origin at the projection center
> Apparently, variant A should be specifiable via a '+no_uoff' flag
> I can specify variant B in this way...
This is a complex topic, see also:
http://trac.osgeo.org/proj/ticket/104
http://trac.osgeo.org/proj/ticket/128
https://trac.osgeo.org/gdal/ticket/4910
On my system running PROJ4 4.9.0 (development version, unreleased), and
using the example from the proj trac ticket 104, I see:
mat <- matrix(c(-(133+(40/60)), 57), ncol=2)
project(mat, '+proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_defs')
[,1] [,2]
[1,] 5e+06 -5e+06
project(mat, '+proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_uoff +no_defs')
[,1] [,2]
[1,] 818676.7 575097.7
but
SPpt <- SpatialPoints(mat, proj4string=CRS("+proj=longlat +datum=WGS84"))
spTransform(SPpt, CRS("+proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_defs"))
SpatialPoints:
coords.x1 coords.x2
[1,] 5e+06 -5e+06
spTransform(SPpt, CRS("+proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_uoff +no_defs"))
SpatialPoints:
coords.x1 coords.x2
[1,] 5e+06 -5e+06
with +no_uoff removed. Indeed:
> CRS("+proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_uoff +no_defs")
CRS arguments:
+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111
+k=0.9999 +x_0=5000000 +y_0=-5000000 +gamma=323.1301023611111
+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
removes +no_uoff.
CRS sends the string out to PROJ4 for checking, via:
> checkCRSArgs("+proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_uoff +no_defs")
[[1]]
[1] TRUE
[[2]]
[1] "+proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_defs"
which drops +no_uoff. The function is:
checkCRSArgs <- function(uprojargs) {
res <- .Call("checkCRSArgs", uprojargs, PACKAGE="rgdal")
res[[2]] <- sub("^\\s+", "", res[[2]])
res
}
and calls the C function:
pj = pj_init_plus(CHAR(STRING_ELT(args, 0)));
then
strcpy(cbuf, pj_get_def(pj, 0));
to get the definition back out. I'm not sure whether pj_get_def is not
dropping this on the way out.
Even subverting the CRS object:
myCRS <- CRS("+proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_uoff +no_defs")
myCRS at projargs <- "+proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_uoff +no_defs"
> spTransform(SPpt, myCRS)
SpatialPoints:
coords.x1 coords.x2
[1,] 5e+06 -5e+06
Coordinate Reference System (CRS) arguments: +proj=omerc +lat_0=57
+lonc=-133.6666666666667 +alpha=323.1301023611111 +k=0.9999
+x_0=5000000 +y_0=-5000000 +gamma=323.1301023611111 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_uoff +no_defs
doesn't get the correct value.
Using proj and cs2cs at the command line give the correct results:
$ proj +proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_defs
133d40'W 57N
5000000.00 -5000000.00
$ proj +proj=omerc +lat_0=57 +lonc=-133.6666666666667
+alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000
+gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_off +no_defs
133d40'W 57N
818676.73 575097.69
$ cs2cs +proj=longlat +datum=WGS84 +to +proj=omerc +lat_0=57
+lonc=-133.6666666666667 +alpha=323.1301023611111 +k=0.9999 +x_0=5000000
+y_0=-5000000 +gamma=323.1301023611111 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs
133d40'W 57N
5000000.00 -5000000.00 0.00
$ cs2cs +proj=longlat +datum=WGS84 +to +proj=omerc +lat_0=57
+lonc=-133.6666666666667 +alpha=323.1301023611111 +k=0.9999 +x_0=5000000
+y_0=-5000000 +gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
+units=m +no_uoff +no_defs
133d40'W 57N
818676.73 575097.69 0.00
Does anyone else have any insights?
Roger
>
> #create a test point, and convert to SPDF
> xy<-c(775082.505359,720337.454359)
> test.pt<-data.frame(x=xy[1],y=xy[2],id=1)
> coordinates(test.pt)<-c('x','y')
>
> #Assign projection
> proj4string(test.pt)<- '+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111
> +k=0.9999 +x_0=5000000 +y_0=-5000000 +gamma=323.1301023611111 +datum=NAD83
> +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0'
>
> # the above is the same as using proj4string(test.pt)<-'+init = epsg:26931',
> # it is also identical to the proj4string obtained for a shapefile created in ArcGIS and imported to R
> # using 'readOGR', that is originally projected in the 26931
>
> #But when I try to specify variant A, the +no_uoff flag is ignored...
> proj4string(test.pt)<- '+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111 +no_uoff
> +k=0.9999 +x_0=5000000 +y_0=-5000000 +gamma=323.1301023611111 +datum=NAD83
> +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0'
> proj4string(test.pt)
> #just to verify, I also try using spTransform...
> test.pt<-spTransform(test.pt,CRS('+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111 +no_uoff
> +k=0.9999 +x_0=5000000 +y_0=-5000000 +gamma=323.1301023611111 +datum=NAD83
> +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0'))
> proj4string(test.pt)
>
> I would report to the os.geo list on this, but I only know the R syntax to this point, so it is difficult to formulate an example for them.
>
> Thanks kindly for any advice you can lend, I appreciate it.
>
> Dave.
>
> ________________________________________________
> sessionInfo()
> R version 3.0.3 (2014-03-06)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> locale:
> [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 LC_NUMERIC=C LC_TIME=English_United States.1252
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] rgdal_0.8-16 raster_2.2-31 sp_1.0-15
>
> loaded via a namespace (and not attached):
> [1] grid_3.0.3 lattice_0.20-27 tools_3.0.3
>
>> library(rgdal)
> rgdal: version: 0.8-16, (SVN revision 498)
> Geospatial Data Abstraction Library extensions to R successfully loaded
> Loaded GDAL runtime: GDAL 1.11.0, released 2014/04/16
> Path to GDAL shared files: C:/Program Files/R/R-3.0.3/library/rgdal/gdal
> GDAL does not use iconv for recoding strings.
> Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
> Path to PROJ.4 shared files: C:/Program Files/R/R-3.0.3/library/rgdal/proj
> _____________________________________________
>
> __________________________________
> Dave Gregovich
> Research Analyst
> Alaska Department of Fish and Game
> Wildlife Conservation Division
> Douglas, AK 99821
> (907) 465-4291
> dave.gregovich at alaska.gov
> __________________________________
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no
More information about the R-sig-Geo
mailing list