| shapefile not in lat/lon [message #86696] |
Thu, 28 November 2013 07:59  |
natha
Messages: 482
Registered: October 2007
|
Senior Member |
|
|
Hi gurus,
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]
Thank you very much in advance,
nata
|
|
|
|
| Re: shapefile not in lat/lon [message #86697 is a reply to message #86696] |
Thu, 28 November 2013 08:22  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
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]
I would set it up like this:
map = cgMap('UTM', Ellipsoid='GRS 1980', CENTER_LONGITUDE=-73.5, $
CENTER_LATITUDE=0, FALSE_EASTING=304800.0)
lonlat = map -> Inverse(xx, yy)
lon = Reform(lonlat[0,*])
lat = Reform(lonlat[1,*])
If you use the Shapefile tools in the Coyote Library it won't matter
much whether your shapefile uses lat/lon values or projected meter
values, as this one does.
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.")
|
|
|
|
|
|
| Re: shapefile not in lat/lon [message #86700 is a reply to message #86698] |
Thu, 28 November 2013 08:50 |
natha
Messages: 482
Registered: October 2007
|
Senior Member |
|
|
No, it is not working or I don't know how to do it... The only thing I want is my shapefile over a map.
This is what I am doing:
cgdisplay, 800, 700, wid=1
map = cgMap('UTM', Ellipsoid='GRS 1980', CENTER_LONGITUDE=-73.5, $
CENTER_LATITUDE=0, FALSE_EASTING=304800.0, LIMIT=[44,-75.25,47,-72.5])
cgmap_grid, /box_axes, map_structure=map
cgmap_continents, /continents, /coasts, /usa, /rivers, /countries, /hires, map_str=map
cgdrawshapes, shapefile, mapcoord=map, /projected_xy
|
|
|
|
| Re: shapefile not in lat/lon [message #86702 is a reply to message #86700] |
Thu, 28 November 2013 08:57 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
nata writes:
>
> No, it is not working or I don't know how to do it... The only thing I want is my shapefile over a map.
> This is what I am doing:
>
> cgdisplay, 800, 700, wid=1
> map = cgMap('UTM', Ellipsoid='GRS 1980', CENTER_LONGITUDE=-73.5, $
> CENTER_LATITUDE=0, FALSE_EASTING=304800.0, LIMIT=[44,-75.25,47,-72.5])
> cgmap_grid, /box_axes, map_structure=map
> cgmap_continents, /continents, /coasts, /usa, /rivers, /countries, /hires, map_str=map
> cgdrawshapes, shapefile, mapcoord=map, /projected_xy
I'm just going out to gather some company for Thanksgiving, but package
all the pieces of the shapefile up in a zip file and send it to me. I'll
have a look when I get the chance. (Good excape from the family!)
But, instead of doing this:
> cgdrawshapes, shapefile, mapcoord=map, /projected_xy
I would try this:
map -> SetProperty, Position=[0.1, 0.1, 0.9, 0.9]
map -> Draw
cgmap_continents, /continents, /coasts, /usa, /rivers, $
/countries, /hires, map_str=map
cgmap_grid, /box_axes, map_structure=map
cgDrawshapes, shapefile
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.")
|
|
|
|
| Re: shapefile not in lat/lon [message #86716 is a reply to message #86696] |
Fri, 29 November 2013 08:17 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
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.")
|
|
|
|
| Transparent Polygon Shapes Displayed on a Map [message #86750 is a reply to message #86716] |
Sun, 01 December 2013 11:18 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> 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. :-)
I learned a lot (and had to update some software!) from this experience.
You can find an article about how I went about displaying a transparent
polygon shape from a shapefile on a map in this article:
http://www.idlcoyote.com/map_tips/transpoly.php
I spent about 15 hours of my time over the past four days figuring this
out and writing an article about it. No one pays me to do this. I do it
because I like to. In the last 10 years I have probably put in 20 hours
a week supporting IDL users here, through e-mail, and on my web page at
no charge to anyone. But, this kind of support for IDL users will be
coming to an end by the end of this year unless I can find some way to
keep doing it.
If this matters to you, if you use the Coyote Library, if you appreciate
the articles on this web page, if I have answered your questions though
e-mail, then I encourage you to take a moment and show your support by
making a donation to this work. I am coming to the reluctant conclusion
that what I do with IDL has no economic value. If this is true, then I
have to find something else to do, and quickly.
A donation in whatever amount you can afford and consider reasonable
will be gratefully appreciated. I love to do this work. It really is the
only thing I have ever wanted to do. I hope it is valuable to you, too.
You can make a donation here if you wish to do so:
http://www.idlcoyote.com/coyotestore/index.php?main_page=ind ex&cPath=67
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
