comp.lang.idl-pvwave archive
Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado

Home » Public Forums » archive » shifted .kml created with IDL
Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend 
Return to the default flat view Create a new topic Submit Reply
Re: shifted .kml created with IDL [message #81030 is a reply to message #80955] Thu, 02 August 2012 00:00 Go to previous messageGo to previous message
titan is currently offline  titan
Messages: 59
Registered: March 2006
Member
On Tuesday, July 31, 2012 6:36:37 PM UTC+2, David Fanning wrote:
> David Fanning writes:
>
>
>
>> I'm not sure if these are the *centers* of the corner
>
>> pixels or the outside edge. I would probably get this
>
>> information directly from the GeoTiff file. If you did,
>
>> you would be getting the outside edge boundaries.
>
>>
>
>> In other words, you can get this information from the
>
>> GeoTiff structure using tie points, the scale in the
>
>> X and Y direction, and the size of the image. All info
>
>> that is available for the GeoTiff file:
>
>>
>
>> http://www.idlcoyote.com/map_tips/pixel_to_ll.html
>
>
>
> All this boundary information is more easily obtained
>
> with the Coyote Library routine FindMapBoundary, of course:
>
>
>
> void = FindMapBoundary(geoTiffFile, b, XRANGE=xr, YRANGE=yr)
>
>
>
> The boundary is in the variable b. Or you can use the ranges.
>
> Same numbers, just presented differently, depending on what
>
> you need.
>
>
>
> Cheers,
>
>
>
> David
>
>
>
>
>
> --
>
> David Fanning, Ph.D.
>
> Fanning Software Consulting, Inc.
>
> Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
>
> Sepore ma de ni thui. ("Perhaps thou speakest truth.")

Hi David,
I tried to collect all your useful suggestions to get the rid of this problem but I think I'm still missing something and I'm getting confused :( (Do I need holiday???)

I had two approaches:
1. I'm calculating the center of the corner pixels and from this (using http://www.idlcoyote.com/map_tips/pixel_to_ll.html) I calculated the lat and lon for all pixel of the image in the Google Earth projection


tiff_image=path_fname+tiff_image

;mapCoord = cgGeoMap(tiff_image)

;geoTiffMapStruct = mapCoord -> GetMapStruct()
geMapStruct = Map_Proj_Init('Mercator', /GCTP, Datum=8)
image=Read_Tiff(tiff_image, GEOTIFF=geotag)

image = Reverse(image, 2) ; Necessary for TIFF files.
s=Size(image, /Dimensions)
xscale = geotag.ModelPixelScaleTag[0]
yscale = geotag.ModelPixelScaleTag[1]
tp = geotag.ModelTiePointTag[3:4]

xOrigin = tp[0]
yOrigin = tp[1] - (yscale * s[1])
xEnd = xOrigin + (xscale * s[0])
yEnd = tp[1]
xCenter = (xEnd - xOrigin) / 2.0 + xOrigin
yCenter = (yEnd - yOrigin) / 2.0 + yOrigin


uvec = Findgen(s[0]) * xscale + xOrigin
vvec = Findgen(s[1]) * yscale + yOrigin

uarray = Rebin(uvec, s[0], s[1])
varray = Rebin(Reform(vvec, 1, s[1]), s[0], s[1])

lonlat_ge= MAP_PROJ_INVERSE(uarray, varray, MAP_STRUCTURE=geMapStruct)
minLat_ge = Min(lonlat_ge[1,*], MAX=maxLat_ge)
minLon_ge = Min(lonlat_ge[0,*], MAX=maxLon_ge)
limit_ge = [minLat_ge, minLon_ge, maxLat_ge, maxLon_ge]

print, limit_ge
39.646844 6.1113937 42.130700 6.1182209

2. Then I tried the FindMapBoundary function to get the border of my tiff image and then using in a correct way the map_proj_image function. After reading the tiff file and cr4eating a warped image using the Google Earth map structure I convert the uvrange in lat/lon to be passed to my kml.

void = FindMapBoundary(tiff_image, b, XRANGE=xr, YRANGE=yr)
print, 'boundary ', b
print, 'XRANGE ',xr
print, 'YRANGE ',yr

tiffCoord = cgGeoMap(tiff_image)

tiffCoord -> SetProperty, XRANGE=xr, YRANGE=yr

MercatorCoord = Obj_New('cgMap', "MERCATOR", DATUM=8)

warped_r = Map_Proj_Image(REFORM(image[0,*,*]), b, $
IMAGE_STRUCTURE=tiffCoord->GetMapStruct(), MAP_STRUCTURE=MercatorCoord->GetMapStruct(), $
UVRANGE=uvout)

warped_g = Map_Proj_Image(REFORM(image[1,*,*]), b, $
IMAGE_STRUCTURE=tiffCoord->GetMapStruct(), MAP_STRUCTURE=MercatorCoord->GetMapStruct(), $
UVRANGE=uvout)


warped_b = Map_Proj_Image(REFORM(image[2,*,*]), b, $
IMAGE_STRUCTURE=tiffCoord->GetMapStruct(), MAP_STRUCTURE=MercatorCoord->GetMapStruct(), $
UVRANGE=uvout)
print,'uvout'
print,uvout

s=SIZE(warped_r, /DIMENSIONS)
warped = BytArr(3, s[0], s[1])
warped[0, *, *] = ROTATE(warped_r,7)
warped[1, *, *] = ROTATE(warped_g,7)
warped[2, *, *] = ROTATE(warped_b,7)


WRITE_PNG, path_fname+'tiff_warped2mercator.png', warped

lonlat_kml= MAP_PROJ_INVERSE(uvout, MAP_STRUCTURE=geMapStruct)

boundary 680317.23 4878686.4 1045117.2 5152286.4
XRANGE 680317.23 1045117.2
YRANGE 4878686.4 5152286.4

uvout
1252403.4 5414120.5 1789982.0 5829831.8

lonlat_kml
11.250532 43.859578
16.079682 46.500001


The problem is that I expected coordinate near to these
ul_lat = 47.653348
ul_lon = 8.9734384
lr_lat = 45.292876
lr_lon = 13.830713
which are the corners of the area in which I'm interest.
Am I still doing everything wrong?? :(

Thanks for your support!!
[Message index]
 
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Previous Topic: Re: max_value/min_value in the x-axis?
Next Topic: IDL User Group Meeting 2012

-=] Back to Top [=-
[ Syndicate this forum (XML) ] [ RSS ] [ PDF ]

Current Time: Wed Oct 15 16:18:58 PDT 2025

Total time taken to generate the page: 1.35866 seconds