[R-sig-Geo] converting lat/long to North American lambert x/y

[Tomislav Hengl]

You should use the proj4 functionality implemented in the rgdal package.
This is what I get:

> library(rgdal)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.7.0dev, released 2008/11/26
Path to GDAL shared files: C:/PROGRA~1/R/R-210~1.0/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.6.1, 21 August 2008
Path to PROJ.4 shared files: C:/PROGRA~1/R/R-210~1.0/library/rgdal/proj
# pan-american projection:
> AEA <- "+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"
> pnts <- data.frame(lat=c(39.953438, 45.271987), long=c(-119.511746,
-97.339382))
> coordinates(pnts) <- ~long+lat
> proj4string(pnts) <- CRS("+proj=longlat +ellps=WGS84")
> spTransform(pnts, CRS(AEA))
SpatialPoints:
           long     lat
[1,] -1970668.1 2126789
[2,]  -105037.3 2476717
Coordinate Reference System (CRS)
arguments: +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
+towgs84=0,0,0

# your projection:
> LCC <- "+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0
+y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
> spTransform(pnts, CRS(LCC))
SpatialPoints:
            long      lat
[1,] -1864716.27 248057.1
[2,]   -98990.79 551251.4
Coordinate Reference System (CRS)
arguments: +proj=lcc +lat_1=20
+lat_2=60 +lat_0=40 +lon_0=-96
+x_0=0 +y_0=0 +ellps=GRS80
+datum=NAD83 +units=m +no_defs
+towgs84=0,0,0


Hmmm, if you set-up the datum to NAD83, then (I think) you should not
get +towgs84=0,0,0 but it looks like these two things are compatible.

HTH

T. Hengl
More examples with US data: http://spatial-analyst.net/book/NGS.R

> Hello List,
>
> Thanks for any help in advance!  I am new to spatial analysis in R and in
> general and so apologize if I am missing something simple.
>
> I have a set of lat/long coordinates that span North America and I would
> like to convert them to utm's (North_America_Lambert_Conformal_Conic).  I
> have attempted to do this using the mapproject command in the mapproj
> library but I keep getting strangely low utm's.  My attempts are
> replicated
> below
>
> My map attributes:
>
>> summary(NAm)
> Object of class SpatialPolygonsDataFrame
> Coordinates:
>         min     max
> r1 -4907279 4031980
> r2 -3453182 5806061
> Is projected: TRUE
> proj4string :
> [+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0
> +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0]
> Data attributes:
> North America
>             1
>
> sample lat/longs:
>
> x <-    c(39.953438, 45.271987, 43.179207, 49.8703, 41.864192, 47.0049,
> 54.727594, 48.719167, 41.338611, 49, 51.363957, 45.187723, 50.203889,
> 50.094167, 27.479949, 51.384852, 58.943633, 52.4, 55.87)
> y <-    c(-119.511746, -97.339382, -112.280316, -88.939133, -121.144867,
> -99.435282, -107.509804, -88.458611, -112.858611, -94.666667, -105.239754,
> -96.130199, -88.758333, -88.337222, -97.314835,    -96.550598, -113.24913,
> -123.05, -115.3275)
>
> coords<-as.data.frame(cbind(x,y))
>
> My code and results:
>
>> x <-    c(39.953438, 45.271987, 43.179207, 49.8703, 41.864192, 47.0049,
> 54.727594, 48.719167, 41.338611, 49, 51.363957, 45.187723, 50.203889,
> 50.094167, 27.479949, 51.384852, 58.943633, 52.4, 55.87)
>> y <-    c(-119.511746, -97.339382, -112.280316, -88.939133, -121.144867,
> -99.435282, -107.509804, -88.458611, -112.858611, -94.666667, -105.239754,
> -96.130199, -88.758333, -88.337222, -97.314835,    -96.550598, -113.24913,
> -123.05, -115.3275)
>>
>> coords<-as.data.frame(cbind(x,y))
>>
>>
>>
>> mapproject((coords), projection = "lambert", param = c(20, 60),
> orientation = c(40,-96,0))
> $x
>  [1] -0.34331799 -0.11436538 -0.27206521  0.01792251 -0.33693782
> -0.13471935
> -0.16832946  0.02711496 -0.28870789 -0.07032821 -0.17079579
> [12] -0.09824454  0.02089270  0.02835444 -0.15848767 -0.08855674
> -0.17102350
> -0.25815791 -0.20128084
>
> $y
>  [1] -1.401480 -1.797972 -1.514287 -1.986024 -1.369876 -1.750728 -1.574392
> -1.997333 -1.509375 -1.851279 -1.624597 -1.823937 -1.990606 -2.001015
> [15] -1.834957 -1.806036 -1.463011 -1.315917 -1.432837
>
> $range
> [1] -0.34331799  0.02835444 -2.00101509 -1.31591678
>
> $error
> [1] 0
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>