Finding pixel values of GeoTIFF image based on lat/lon (ENVI and IDL give different results).

[Sobriquet]

Hi,

I know this is a relatively simple task but I'm still new to IDL and have been unable to find an appropriate answer to the problem.

I have ~700 GeoTIFF images of dimensions around 4000x4000.I have the lat/long coordinates of 9 points from which I need to extract information. The images are projected on the Albers Conical Equal Area projection. I need to find the closest pixel value (sample, line image coordinates) that matches the lat/lon pairs. I previously used the approach described here: http://www.idlcoyote.com/map_tips/pixel_to_ll.html

However, I was unable to find a matching pair of lat/long in both arrays or convert the closest match between the two to pixel coordinates.
I, then, tried searching through the images using this code:

RES = 100./1.

IMGDATA = READ_TIFF(data_FILES[I], IMAGE_INDEX = 0, GEOTIFF = geotag)
S = SIZE(IMGDATA, /DIMENSIONS); image size

FILE_POS_data = STRSPLIT(data_FILES[I], '\', /EXTRACT, COUNT = POS_data)
FILE_NAME = STRMID(FILE_POS_data[POS_data-1], 0, 15)

RESOLUTION = geotag.MODELPIXELSCALETAG[0] / RES
xmin = geotag.MODELTIEPOINTTAG[3] / RES
ymax = geotag.MODELTIEPOINTTAG[4] / RES
ymin = ymax - (N_ELEMENTS(IMGDATA[0,0,*])) * resolution

PTLAT = my_lat
PTLON = my_lon

alberMap = MAP_PROJ_INIT('Albers Equal Area Conic', $
DATUM=8, $ ; WGS84
SPHERE_RADIUS = 6378137.0, $ ; Earth's equatorial radius, or semi-major axis according to WGS-84
CENTER_LAT=geotag.PROJNATORIGINLATGEOKEY, $
CENTER_LON=geotag.PROJNATORIGINLONGGEOKEY, $
STANDARD_PAR1=geotag.PROJSTDPARALLEL1GEOKEY, $
STANDARD_PAR2=geotag.PROJSTDPARALLEL2GEOKEY)

result = MAP_PROJ_FORWARD(ptlon, ptlat, MAP_STRUCTURE=alberMap)
ptx = result[0, *]/RES
pty = result[1, *]/RES

x = FLOOR((ptx - xmin ) / resolution)
y = -(FLOOR((pty - ymin) / resolution))

PRINT, x,y,resolution

results = IMGDATA[x-1,y-1]

The problem is that when I double check the results using the same latitude and longitude in ENVI I get different results/pixels.

For example, when I input the my_lat/my_lon pair (64.7000, -148.23300) in the IDL code, I get a pixel value of 3164 and a sample/line value of 805, 2501. However, when I open the image in ENVI and change the projection to lat/lon, I get a pixel value of 3122, and a sample/line value of 813,2495. This might be different because ENVI could be using a different reprojection scheme than IDL.

When I check for the lat/lons that get approximated for the sample/line combinations in ENVI, I get that:

My_lat, mylon = 64.7000, -148.23300

ENVI 813, 2495 = 64.7002, -148.2349
IDL 805, 2501 = 64.6955, -148.2527

ENVI gives me more precise results, and although a difference between 3164 and 3122 is relatively small, other images and coordinate pairs give me more wildly varying results (some have a difference of 1500, 2000).

My questions are these:

Is there any method to reconcile this discrepancy in my IDL code?
What is the reprojection scheme that ENVI uses to change from ALBERS to lat/lon, and that could explain the difference?

Any suggestions would be greatly appreciated.

Thanks.

[Coyote]

Post by Sobriquet
I have ~700 GeoTIFF images of dimensions around 4000x4000.I have the lat/long coordinates of 9 points from which I need to extract information. The images are projected on the Albers Conical Equal Area projection. I need to find the closest pixel value (sample, line image coordinates) that matches the lat/lon pairs. I previously used the approach described here: http://www.idlcoyote.com/map_tips/pixel_to_ll.html
However, I was unable to find a matching pair of lat/long in both arrays or convert the closest match between the two to pixel coordinates.

Humm. I'm not sure you are following the approach described in that article, because I see no evidence in your code that you are reversing your TIFF image in Y, etc., etc.

If I were going to do this I *would* follow the approach in the article, exactly, up to the point where I had the two vectors uvec and vvec. Then,
I would convert the point (presumably in lat/lon) you want to find into projected meter values and find the image index with Value_Locate, like this:

pt_lon = -148.23300
pt_lat = 64.7000

xy = Map_Proj_Forward(pt_lon, pt_lat, Map_Structure=Albermap)
pt_x = xy[0]
pt_y = xy[1]

xindex = Value_Locate(uvec, pt_x)
yindex = Value_Locate(vvec, pt_y)

Print, 'Image Value: ', image[xindex,yindex]

ll = Map_Proj_Inverse(uvec[xindex], vvec[yindex], Map_Structure=Albermap)
Print, 'Nearest Pixel Location (lon/lat): '
Print, ' Longitude: ', ll[0], ' Image X Coord: ', xindex
Print, ' Latitude: ', ll[1], ' Image Y Coord: ', yindex

I wouldn't use the SPHERE_RADIUS keyword in Map_Proj_Init either. I can't tell if it is hurting you, but it can't be doing you any good. :-)

Cheers,

David

[Sobriquet]

David,

Thanks for the reply.

Your suggestion to remove the SPHERE_RADIUS keyword was accurate: removing it makes no difference.

Here is what I got:

Image Value: 3085.77
Nearest Pixel Location (lon/lat):
Longitude: -148.23405 Image X Coord: 805
Latitude: 64.699193 Image Y Coord: 1662

What's interesting is that when I input the lat/lons in ENVI it approximates it to the same values that I get when I enter my original lat/lons (64.7000, -148.23300 // 3122). However, when I input the x,y coordinates into ENVI, I get this

(65.4425, -148.09711726 // 3552.166016),

which is very far from my original location. In fact,when I enter the nearest pixel location lat/lons into Google Earth, it lands ~85km away from my input lat/long.

Now, I know ENVI is not perfect. A quick survey on Google Earth also shows that the approximation of the location I entered is about 100 m off.

I think the problem is that since the resolution of the arrays is very high (100 m), there are a lot of possible matches in the position vectors. Value_Locate is searching through the u and v vectors individually and finding the closest match for each value without taking into account the relationship between the two vectors,i.e., it is not finding the closest u match in relation to the v match that together make the nearest point possible to the input values. As a result, you might get very close u and v individual matches for your input lat/lon that translate into a very off pair of sample/line values in your pixel coordinates.

I hate to be a stickler, but is there a way to establish a relationship between the two position vectors that would help solve the ambiguity?

Thanks again

[David Fanning]

Post by Sobriquet
Image Value: 3085.77
Longitude: -148.23405 Image X Coord: 805
Latitude: 64.699193 Image Y Coord: 1662
What's interesting is that when I input the lat/lons in ENVI it approximates it to the same values that I get when I enter my original lat/lons (64.7000, -148.23300 // 3122). However, when I input the x,y coordinates into ENVI, I get this
(65.4425, -148.09711726 // 3552.166016),
which is very far from my original location. In fact,when I enter the nearest pixel location lat/lons into Google Earth, it lands ~85km away from my input lat/long.

Well, ENVI, remember, uses a different coordinate system than
IDL does. Theirs starts at 1, IDL's starts at 0.

Post by Sobriquet
Now, I know ENVI is not perfect. A quick survey on Google Earth also shows that the approximation of the location I entered is about 100 m off.

Did you convert your coordinate from an Albers projection to the
Equirectangular projection used by Google Earth? If not, this
can account for differences of at least 100 meters, depending
upon where you are on the Earth.

Post by Sobriquet
I think the problem is that since the resolution of the arrays is very high (100 m), there are a lot of possible matches in the position vectors. Value_Locate is searching through the u and v vectors individually and finding the closest match for each value without taking into account the relationship between the two vectors,i.e., it is not finding the closest u match in relation to the v match that together make the nearest point possible to the input values. As a

result, you might get very close u and v individual matches for your input lat/lon that translate into a very off pair of sample/line values in your pixel coordinates.

On the contrary, there is only one possible match on this grid. It is
true that the actual location of the grid crossing point can be off from
the true location of your point by at something like 100 meters (if this
is the size of the grid), but there are not "multiple" matches. There is
one possible match to a particular point.

Post by Sobriquet
I hate to be a stickler, but is there a way to establish a relationship between the two position vectors that would help solve the ambiguity?

Value_Locate is finding the index of the value that is closest to your
chosen point that doesn't exceed your point value. You could subtract
the grid value from your point, and if it is larger than half a grid
size, you could say the next index is closer than the one Value_Locate
found.

xindex = Value_Locate(uvec, pt_x)
IF pt_x = uvec[xindex] > (uvec[1]-uvec[0])/2 THEN xindex = xindex + 1

I guess it depends on how anal you want to be. :-)

Cheers,

David

[David Fanning]

Post by Coyote
xindex = Value_Locate(uvec, pt_x)
IF pt_x = uvec[xindex] > (uvec[1]-uvec[0])/2 THEN xindex = xindex + 1

Whoops!

Shouldn't write code with one hand on a burrito! :-(

xindex = Value_Locate(uvec, pt_x)
IF (pt_x - uvec[xindex]) GT ((uvec[1]-uvec[0])/2.0) THEN $
xindex = xindex + 1

Cheers,

David

[Sobriquet]

Thanks for the help!

I found that indeed the subtraction between the grid value from my point was more than 50 (~90). I'll go with the next pixel over.