Skip to content

function to convert to/from U.S. State Plane Coordinates?

15 messages · Waichler, Scott R, Hodgess, Erin, Rob Robinson +2 more

Original
Waichler, Scott R
# 
Hi,  Does an R function exist to convert to and from the State Plane Coordinates system used in the U.S.?  My search with Rseek was fruitless.

Scott Waichler
Hydrology Group, Energy & Environment Directorate
Pacific Northwest National Laboratory
Richland, WA, USA
scott.waichler at pnnl.gov
Hodgess, Erin
# 
Could you mean spTransform from rgdal, please?
Waichler, Scott R
# 
Erin, thanks, I'm getting closer.  Rgdal's projInfo() lists UTM, which I believe is closely related to the State Plane Coordinate system, but I still don't see a direct transformation to SPC in R.  --Scott Waichler
Michael Sumner
# 
State plane is either UTM or Lambert Conformal Conic ("+proj=lcc"),
but this is a *family* (of families) of projections, not a specific
coordinate system. You'll need to be more specific by finding from an
authority which actual "State Plane" projection you need.

http://en.wikipedia.org/wiki/State_Plane_Coordinate_System

spatialreference.org is a good place to look if you are lacking detail:

http://spatialreference.org/ref/?search=state+plane

Often there is metadata present in data files so if you have those it
might be worth describing them, and trying tools in rgdal to find this
out. An EPSG code likely to be the most concise identifier, but it
could be PROJ.4 or WKT or something else that your colleagues might
know.

Cheers, Mike.


On Tue, Sep 17, 2013 at 9:09 AM, Waichler, Scott R
<Scott.Waichler at pnnl.gov> wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Waichler, Scott R
    

    

      
      #
    

  


  

      
Ok, now I know about the existence of EPSG codes, but apparently my rgdal installation can't find the list of them?  Do I need to find a certain file in my system's PROJ4 installation, then tell that to rgdal?

            

      
      
      
    

        
Error in make_EPSG() : Error opening epsg file

            

      
      
      
    

        
Error in spTransform(pnw, CRS = CRS("+init=epsg:102749")) :
  error in evaluating the argument 'CRSobj' in selecting a method for function 'spTransform': Error in CRS("+init=epsg:102749") : no system list, errno: 13

--Scott

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Hodgess, Erin
    

    

      
      #
    

  


  

      
Hi Scott:

Are you in Windows or a Linux type computer, please?



-----Original Message-----
From: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Waichler, Scott R
Sent: Monday, September 16, 2013 6:49 PM
To: Rob Robinson
Cc: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] function to convert to/from U.S. State Plane Coordinates?

Ok, now I know about the existence of EPSG codes, but apparently my rgdal installation can't find the list of them?  Do I need to find a certain file in my system's PROJ4 installation, then tell that to rgdal?

            

      
      
      
    

        
Error in make_EPSG() : Error opening epsg file

            

      
      
      
    

        
Error in spTransform(pnw, CRS = CRS("+init=epsg:102749")) :
  error in evaluating the argument 'CRSobj' in selecting a method for function 'spTransform': Error in CRS("+init=epsg:102749") : no system list, errno: 13

--Scott

            

      
      
      
    

        
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

  
  


  


      

    

    

      
      

        


  

        

      Waichler, Scott R
    

    

      
      #
    

  


  

      
I am working in Linux.  

I clicked on the ESRI:102749 option on the spatialreference.org website that Rob referred to, and got a spec that looks like it should work inside CRS(), but something seems to be amiss.

            

      
      
      
    

        
Error in spTransform(pnw, CRS = CRS("+proj=lcc +lat_1=45.83333333333334 +lat_2=47.33333333333334 +lat_0=45.33333333333334 +lon_0=-120.5 +x_0=500000.0000000002 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs")) :
  load package rgdal for spTransform methods

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Hodgess, Erin
    

    

      
      #
    

  


  

      
I think that you need to download some proj4 stuff.  If you're on Ubuntu, go to the Ubuntu package search page and you'll find it.

Also, what is the str(pnw.spc) please?

Thanks,
Erin


-----Original Message-----
From: Waichler, Scott R [mailto:Scott.Waichler at pnnl.gov] 
Sent: Monday, September 16, 2013 7:10 PM
To: Hodgess, Erin; Rob Robinson
Cc: r-sig-geo at r-project.org
Subject: RE: [R-sig-Geo] function to convert to/from U.S. State Plane Coordinates?

I am working in Linux.  

I clicked on the ESRI:102749 option on the spatialreference.org website that Rob referred to, and got a spec that looks like it should work inside CRS(), but something seems to be amiss.

            

      
      
      
    

        
Error in spTransform(pnw, CRS = CRS("+proj=lcc +lat_1=45.83333333333334 +lat_2=47.33333333333334 +lat_0=45.33333333333334 +lon_0=-120.5 +x_0=500000.0000000002 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs")) :
  load package rgdal for spTransform methods

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Hodgess, Erin
    

    

      
      #
    

  


  

      
Whoops...I meant pnw


-----Original Message-----
From: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Hodgess, Erin
Sent: Monday, September 16, 2013 7:12 PM
To: Waichler, Scott R; Rob Robinson
Cc: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] function to convert to/from U.S. State Plane Coordinates?

I think that you need to download some proj4 stuff.  If you're on Ubuntu, go to the Ubuntu package search page and you'll find it.

Also, what is the str(pnw.spc) please?

Thanks,
Erin


-----Original Message-----
From: Waichler, Scott R [mailto:Scott.Waichler at pnnl.gov]
Sent: Monday, September 16, 2013 7:10 PM
To: Hodgess, Erin; Rob Robinson
Cc: r-sig-geo at r-project.org
Subject: RE: [R-sig-Geo] function to convert to/from U.S. State Plane Coordinates?

I am working in Linux.  

I clicked on the ESRI:102749 option on the spatialreference.org website that Rob referred to, and got a spec that looks like it should work inside CRS(), but something seems to be amiss.

            

      
      
      
    

        
Error in spTransform(pnw, CRS = CRS("+proj=lcc +lat_1=45.83333333333334 +lat_2=47.33333333333334 +lat_0=45.33333333333334 +lon_0=-120.5 +x_0=500000.0000000002 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs")) :
  load package rgdal for spTransform methods

            

      
      
      
    

        

      
      
    

        
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

  
  


  


      

    

    

      
      

        


  

        

      Waichler, Scott R
    

    

      
      #
    

  


  

      
My system is Redhat, RHEL5.

            

      
      
      
    

        
List of 4
 $ x    : num [1:893] -117 -117 -117 -117 -117 ...
 $ y    : num [1:893] 46 46 45.9 45.9 45.9 ...
 $ range: num [1:4] -125 -111 42 49
 $ names: chr [1:7] "idaho" "oregon" "washington:san juan island" "washington:lopez island" ...
 - attr(*, "class")= chr "map"

So, it sounds like there is extra proj4 stuff that I need to install beyond the base installation I already did, which also enabled installation of rgdal? 

Scott

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Hodgess, Erin
    

    

      
      #
    

  


  

      
Ok.  I think I see it now.

You need to have a spatial object rather than a list.  

So you might set up 

	pnw1 <- data.frame(x=pnw$x,y=pnw$y)
	coordinates(pwn1) <- ~x+y
Then you'll have a spatial object.  After that, put in your projection:
	proj4string(pnw1) <- CRS(xxxxx)
	pnw2 <- spTransform(pnw1,CRS(yyyy))




-----Original Message-----
From: Waichler, Scott R [mailto:Scott.Waichler at pnnl.gov] 
Sent: Monday, September 16, 2013 7:19 PM
To: Hodgess, Erin; Rob Robinson
Cc: r-sig-geo at r-project.org
Subject: RE: [R-sig-Geo] function to convert to/from U.S. State Plane Coordinates?

My system is Redhat, RHEL5.

            

      
      
      
    

        
List of 4
 $ x    : num [1:893] -117 -117 -117 -117 -117 ...
 $ y    : num [1:893] 46 46 45.9 45.9 45.9 ...
 $ range: num [1:4] -125 -111 42 49
 $ names: chr [1:7] "idaho" "oregon" "washington:san juan island" "washington:lopez island" ...
 - attr(*, "class")= chr "map"

So, it sounds like there is extra proj4 stuff that I need to install beyond the base installation I already did, which also enabled installation of rgdal? 

Scott

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Tue, 17 Sep 2013, Waichler, Scott R wrote:

        

            

      
      
      
    

        
The package does document all of these requirements:

file.show(system.file("README", package="rgdal"))

or:

https://r-forge.r-project.org/scm/viewvc.php/pkg/inst/README?view=markup&root=rgdal

or read the inst/README file from the source package tarball. Note that 
the Fedora/RHEL case is covered specifically.

Nothing about projections is ever simple, nobody should assume that it is. 
The assumption is very anachronistic, as clear standards have only been 
available recently (say < 20 years, even < 10 years for GNSS/GPS), and 
most maps are older; mapping agencies value continuity, and use their 
legacy definitions, which vary greatly from place to place, and are often 
not well defined.

Note also the link on the Spatial Task View to:

http://spatialreference.org/

as a resource for investigating which incantation may match your data. As 
has been mentioned, State Plane contains hundreds of candidate 
definitions, depending among other things on the part of the state 
concerned and the datum in use at the time the data were created. 
Searching this website for "Washington" - your location state - gives 18 
hits in EPSG, and only you can know which is appropriate.

Hope this clarifies,

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Waichler, Scott R
    

    

      
      #
    

  


  

      
Erin,

Thanks, now I see that I was trying to use spTransform() on a map object instead of a Spatial one.  I followed your logic and came up with the following to plot the state of Washington using state plane coordinates instead of the original geographic values.  But I get an extraneous triangular polygon.  I am guessing it has to do with not repeating the initial x,y at the end of a polygon, or a missing row of NA.  What am I missing? 

library(maps)
library(rgdal)
wa <- map('state', region = 'washington', plot=F, type="n")
ind.not.na <- which(!is.na(wa$x))  # note which rows are not NA
wa1 <- data.frame(x=wa$x,y=wa$y)  # make a Spatial object from the map object
coordinates(wa1) <- ~x+y
proj4string(wa1) <- CRS("+proj=longlat +datum=WGS84")
# I got the CRS spec from http://spatialreference.org/ref/esri/102749/proj4/
wa.spc <- spTransform(wa1, CRS("+proj=lcc +lat_1=45.83333333333334 +lat_2=47.33333333333334 +lat_0=45.33333333333334 +lon_0=-120.5 +x_0=500000.0000000002 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs"))

# make another data frame with the same number of rows as the map() original
wa2 <- data.frame(x=wa$x,y=wa$y)
# repace the geographic values with the state plane coordinate values
wa2$x[ind.not.na] <- wa.spc$x  
wa2$y[ind.not.na] <- wa.spc$y
plot(wa2, type="n")
polygon(wa2)  # gives a spurious triangular polygon--why?

            

      
      
      
    

        
--Scott

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Tue, 17 Sep 2013, Waichler, Scott R wrote:

        

            

      
      
      
    

        
You have not looked for information on converting legacy map objects to sp 
classes. Look for ?map2SpatialPolygons in the maptools package. Better 
still, find a better representation of the state boundary as a shapefile 
and use that instead.

Roger