On Wed, Nov 17, 2010 at 10:44 AM, Caroline Keef <
caroline.k...@jbaconsulting.co.uk> wrote:
> The shapefile was created in ARC and according to the .prj the details of
> the projection are
>
>
>
> PROJCS["North_Pole_Lambert_Azimuthal_Equal_Area",
>
> GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",
>
> SPHEROID["WGS_1984",6378137.0,298.257223563]],
>
> PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
>
> PROJECTION["Lambert_Azimuthal_Equal_Area"],
>
> PARAMETER["False_Easting",0.0],
>
> PARAMETER["False_Northing",0.0],
>
> PARAMETER["Central_Meridian",0.0],
>
> PARAMETER["Latitude_Of_Origin",90.0],
>
> UNIT["Meter",1.0]]
>
>
>
> According to ARC map the coordinates of Mumbai should be (roughly) 7048812,
> -2172723 the answers I get are
>
>
>
> mapproject(x=72.83,y=18.975,proj="azequalarea")
>
> $x
>
> [1] 0
>
>
>
> $y
>
> [1] -0.8214892
>
>
>
> $range
>
> [1] 0.0000000 0.0000000 -0.8214892 -0.8214892
>
>
>
> $error
>
> [1] 0
>
>
>
>
>
> project(xy=c(72.83,18.975),"+proj=laea +lat_0 = -90 +lon_0=0 +x_0=0 +y_0 =
> 0")
>
> [,1] [,2]
>
> [1,] 7197711 2590300
>
> > [,1] [,2]
>
>
>
> > project(xy=c(72.83,18.975),"+proj=laea +lat_0 = 0 +lon_0=0 +x_0=0 +y_0 =
> 0")
>
> [,1] [,2]
>
> [1,] 7197711 2590300
>
>
I wouldn't use mapproject - get sp and rgdal and use the spTransform
function. Like this:
d=data.frame(x=72.83,y=18.975,p="mumbai")
coordinates(d)=~x+y
- now I have a spatial points data frame. This next bit might be wrong:
proj4string(d)="+init=epsg:4326"
- since you didn't say what coordinate system the lat-long for Mumbai was.
epsg:4326 is WGS84 (unprojected).
Anyway, let's transformobulate:
spTransform(d,CRS("+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0"))
coordinates p
1 (7078530, -2187110) mumbai
- this is close to the numbers you gave (these units are metres, so decide
for yourself what 'close' means! A negative "northing" is because your false
origin is at the north pole (+90). I think George Bush once referred to
India as being in the southern hemisphere...
I note you had +lat_0=-90 in your code but your text gave +90. Anyways,
spTransform does give different answers:
> spTransform(d,CRS("+proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0"))
coordinates p
1 (9904750, 3060350) mumbai
> spTransform(d,CRS("+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0"))
coordinates p
1 (7078530, -2187110) mumbai
> spTransform(d,CRS("+proj=laea +lat_0=0 +lon_0=0 +x_0=0 +y_0=0"))
coordinates p
1 (7208810, 2576930) mumbai
so it's doing _something_.
So, it might be necessary to do some more tuning of the CRS() projection
parameters, otherwise, that's how I'd do it!
Barry
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
