| Re: lon-lat coordinates for Geotiff images? [message #72034] |
Wed, 04 August 2010 05:49  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Dave Poreh writes:
> I have some geotiff images (DEM). I can read these files with:
> im=read_tiff(file)
> And I can contour them:
> CONTOUR, im, /FILL, NLEVELS = 7
> But in the contour I get the pixel locations in the x-y axis (local
> coordinates). Does anybody know how could I get lon-lat coordinates in
> x-y axis?
In a contour plot you have some kind of regular (or irregular)
grid. That is to say, your labels are labels of straight lines.
That is seldom true of a map, where straight lines don't seem
to exist.
There is a grid of straight lines associated with a map-
projected image, however. This is the XY grid, which can
be read directly from the GeoTiff file:
http://www.dfanning.com/map_tips/tiffoverlay.html
I would create these vectors using tools from the
Catalyst Library:
mapCoord = GeoCoord('mygeoTiffFile.tif')
mapCoord -> GetProperty, XRANGE=xr, YRANGE=yr
image = Read_Tiff(('mygeoTiffFile.tif')
s = Size(image, /DIMENSIONS)
xvec = Scale_Vector(Findgen(s[0]), xr[0], xr[1])
yvec = Scale_Vector(Findgen(s[1]), yr[0], yr[1])
Then, you could draw your contour plot:
pos = [0.15, 0.15, 0.9, 0.9]
Contour, Reverse(image,2), xvec, yvec, XStyle=1, $
YStyle=1, Position=pos, NLEVELS=7, /FILL
If you wanted to know the latitude and longitude of
any point in the image, you could convert the XY location
to a lon/lat location:
mapStruct = mapCoord->GetMapStructure()
lonlat = Map_Proj_Inverse(x, y, MAP_STRUCTURE=mapStruct)
lon = lonlat[0]
lat = lonlat]1]
If you wanted to see a lat/lon map grid on your contour map, you
could easily add it like this:
grid = Obj_New('Map_Grid', MAPCOORD=mapCoord, /AUTODRAW_GRID, $
COLOR='gray')
mapCoord -> SetProperty, POSITION=pos, GRID_OBJECT=grid
mapCoord -> Draw, /OVERLAYS
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
| Re: lon-lat coordinates for Geotiff images? [message #72157 is a reply to message #72034] |
Thu, 05 August 2010 02:39  |
d.poreh
Messages: 406
Registered: October 2007
|
Senior Member |
|
|
On Aug 4, 5:49 am, David Fanning <n...@dfanning.com> wrote:
> Dave Poreh writes:
>> I have some geotiff images (DEM). I can read these files with:
>> im=read_tiff(file)
>> And I can contour them:
>> CONTOUR, im, /FILL, NLEVELS = 7
>> But in the contour I get the pixel locations in the x-y axis (local
>> coordinates). Does anybody know how could I get lon-lat coordinates in
>> x-y axis?
>
> In a contour plot you have some kind of regular (or irregular)
> grid. That is to say, your labels are labels of straight lines.
> That is seldom true of a map, where straight lines don't seem
> to exist.
>
> There is a grid of straight lines associated with a map-
> projected image, however. This is the XY grid, which can
> be read directly from the GeoTiff file:
>
> http://www.dfanning.com/map_tips/tiffoverlay.html
>
> I would create these vectors using tools from the
> Catalyst Library:
>
> mapCoord = GeoCoord('mygeoTiffFile.tif')
> mapCoord -> GetProperty, XRANGE=xr, YRANGE=yr
> image = Read_Tiff(('mygeoTiffFile.tif')
> s = Size(image, /DIMENSIONS)
> xvec = Scale_Vector(Findgen(s[0]), xr[0], xr[1])
> yvec = Scale_Vector(Findgen(s[1]), yr[0], yr[1])
>
> Then, you could draw your contour plot:
>
> pos = [0.15, 0.15, 0.9, 0.9]
> Contour, Reverse(image,2), xvec, yvec, XStyle=1, $
> YStyle=1, Position=pos, NLEVELS=7, /FILL
>
> If you wanted to know the latitude and longitude of
> any point in the image, you could convert the XY location
> to a lon/lat location:
>
> mapStruct = mapCoord->GetMapStructure()
> lonlat = Map_Proj_Inverse(x, y, MAP_STRUCTURE=mapStruct)
> lon = lonlat[0]
> lat = lonlat]1]
>
> If you wanted to see a lat/lon map grid on your contour map, you
> could easily add it like this:
>
> grid = Obj_New('Map_Grid', MAPCOORD=mapCoord, /AUTODRAW_GRID, $
> COLOR='gray')
> mapCoord -> SetProperty, POSITION=pos, GRID_OBJECT=grid
> mapCoord -> Draw, /OVERLAYS
>
> Cheers,
>
> David
> --
> David Fanning, Ph.D.
> Fanning Software Consulting, Inc.
> Coyote's Guide to IDL Programming:http://www.dfanning.com/
> Sepore ma de ni thui. ("Perhaps thou speakest truth.")
Thanks David it was very helpful!
Cheers
Dave
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
