Prev: wavelet coefficients for data clustering
Next: supervised learning of semantic classes for image annotation and retrieval source code
From: Rossella on 27 Apr 2010 13:30
Hi,
I am importing some shapefiles into Matlab and I need to retrieve the latidude and longitude coordinates of them (currently in planar coordinates). From the projection file (.prj) I accessed the type of projection (Transverse Mercator) and its parameters (in the following the text of the .prj file).

PROJCS["NAD_1983_StatePlane_Illinois_East_FIPS_1201_Feet",
GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",984250.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",-88.33333333333333],
PARAMETER["Scale_Factor",0.999975],
PARAMETER["Latitude_Of_Origin",36.66666666666666],
UNIT["Foot_US",0.3048006096012192]]

Based on these information I created the projection structure:

mstruct = defaultm('tranmerc');
mstruct.geoid = [6378137.0 0.08181919];
mstruct.origin = [36.667 -88.333 0];
mstruct.falseeasting = 984250.0;
mstruct.falsenorthing = 0.0;
mstruct.scalefactor = 0.999975;
mstruct = defaultm(mstruct);

and then I transformed the coordinates from planar to geographic (lat & lon) with minvtran:

[lat lon] = minvtran(mstruct, ShapeFile.X, ShapeFile.Y);

The problem is that the resulting lat/long coordinates are not correct (the latidudes points are shifted of about 11° North, while the longitude are shifted of about 2° East). I tried to modify some of the parameters I am not sure about (eccentricity and false easting), and I obtained new lat/long coordinates, but still not correct. I tried also the command projinv, but I obtained the same results.
I do not understand where I'm doing wrong.

Any help or suggestion is very appreciated.
Thanks in advance.
Rossella



From: Rob Comer on 4 May 2010 17:13
You're probably having a problem with units. It seems as if everything is in U.S. Survey Feet except for the semimajor axis of the reference ellipsoid (mstruct.geoid), which is in meters. Try converting like this: 6378137 * unitsratio('survey foot','meter') or get the value from directly almanac: almanac('earth','grs80','survey foot').

Rob Comer
Mapping Toolbox Development
MathWorks
From: Rossella on 4 May 2010 17:48
Thank you! I thought I had to put meters as unit because the radius of the earth in the projection file was expressed in meters.

Rossella
 | 
Pages: 1
Prev: wavelet coefficients for data clustering
Next: supervised learning of semantic classes for image annotation and retrieval source code