On Wed, Nov 17, 2010 at 10:44 AM, Caroline Keef <
caroline.keef@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]]

