nata writes:
> I have a set of shapefiles defining some regions but the content is not in lat/lon coordinates. I guess I should use the projection defined in the prj file but I don't exactly know how to do it.
> I am hoping for some help here.
>
> This is the content of the prj file:
> PROJCS["NAD_1983_MTM_8",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",304800.0],PARAMETER["False_Northing ",0.0], PARAMETER["Central_Meridian",-73.5],PARAMETER["Scale_Factor ",0.9999], PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]
>
> Using this information how could I get the lat/lon values of the area defined as:
>
> xx=[314000.00,314000.00,315000.00,315000.00,314000.00]
> yy=[4987000.0,4988000.0,4988000.0,4987000.0,4987000.0]
Nata sent me this file, and I was *finally* able to figure out what is
going on. This file is a Modified UTM projection used in Eastern Canada.
It's zone (zone 8 in this case) is only 3 degrees wide, rather than the
usual 6 degrees wide. Initially, I tried to use a UTM map projection
with the equivalent zone number (18 in this case), but the resulting
shape was just a little bit out of place.
I fooled around with the UTM projection for quite some time, and then I
realized that whatever value I was using for FALSE_EASTING (304800.0 in
this case) made absolutely no difference when I converted the projected
XY values in the shapefile to latitude and longitude. This turns out to
be because the UTM projection is *assuming* the usually 5000000.0 false
easting value for this projection. It ignores whatever you use with the
FALSE_EASTING keyword!
This was the clue I needed to realize I needed to change the map
projection. I chose a Transverse Mercator projection, with the proper
FALSE_EASTING value, and suddenly all appeared well! This also had the
advantage that I could use more of the information from the file. For
example, now the CENTER_LONGITUDE and MERCATOR_SCALE values made sense!
utmmap = cgMap('Transverse Mercator', Ellipsoid='GRS 1980', $
MERCATOR_SCALE=0.9999, FALSE_EASTING=304800., $
CENTER_LONGITUDE=-73.5)
I wanted to put this on a Google map so I could make sure it was located
correctly, so then it was simply a matter of converting these Transverse
Mercator coordinates, using the GRS 1980 datum, into the Google Mercator
coordinates, with the WGS84 datum, of the Google map:
roi = cgExtractShape(shapefile, 'LAYER', 'RICHELIEU')
roi -> GetProperty, DATA=xy
ll = Map_Proj_Inverse(xy[0,*], xy[1,*], MAP=utmmap->GetMapStruct())
xy = Map_Proj_Forward(ll[0,*], ll[1,*], MAP=googlemap->GetMapStruct())
roi -> SetProperty, DATA=xy
cgDraw_ROI, roi, /Outline, Color='red', Thick=4
Wahla! Perfect. :-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thue. ("Perhaps thou speakest truth.")
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
