[R-sig-Geo] a projection question

[Hodgess, Erin]

Great!  This is so cool!

Thanks!

________________________________________
From: b.rowlingson at gmail.com [b.rowlingson at gmail.com] on behalf of Barry Rowlingson [b.rowlingson at lancaster.ac.uk]
Sent: Sunday, December 22, 2013 4:16 AM
To: Hodgess, Erin
Cc: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] a projection question

First thought is to try and project a sample point and catch any error
messages, but you can't be sure that your sample point is valid in the
given coordinate system.

Second thought was to look at the proj4 and rgdal packages to see if
they implemented checking of PROJ.4 strings. Couldn't see it.

So third thought was to look at the PROJ.4 API and see what its got:

http://trac.osgeo.org/proj/wiki/pj_init_plus

which is a C function to initialise a projection. So I wrote a little wrapper:

#include <proj_api.h>

void ptest(char *proj[], int *status){
  projPJ out;
  if(!(out = pj_init_plus(proj[0]))){
    *status = 1;
  }else{
    *status = 0;
   /* need to pj_free the out object here... I think... */
  }
}

with an R wrapper: isValidProjection <-
function(s){.C("ptest",s,status=as.integer(-1))$status==0}

compiled it, linked with libproj.so and tested it:


> isValidProjection("+init=epsg:4326")
[1] TRUE

- user typo:

> isValidProjection("+init=apsg:4326")
[1] FALSE

- try more than +init strings:

> isValidProjection("+proj=tmerc +lon_0 +datum=WGS84")
[1] TRUE

- some bogus EPSG codes:

> isValidProjection("+init=epsg:4326999")
[1] FALSE
> isValidProjection("+init=epsg:-1")
[1] FALSE

Looks like it works. The tricky part of all this is linking the C code
with libproj.so which I did by doing:

R CMD SHLIB projtest.c /usr/lib/libproj.so

on the command line. Might work for you, but it might be worth asking
Roger, as maintainer or rgdal, or Smon Urbanek as maintainer of
project, to add it as a top-level R function.

 Maybe I'll check out the rgdal source and push a change for Roger...

Barry




On Sat, Dec 21, 2013 at 10:19 PM, Hodgess, Erin <HodgessE at uhd.edu> wrote:
> Hello everyone:
>
> I'm setting up a function in which the user enters the projection (as part of other stuff).
>
> What is the best way to determine if the entered projection is correct, please? (so I can send back an error message right away).
>
> Thanks,
> Erin
>
>
>         [[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