Folks,
OK, I have warped my GOES image into an Albers Equal Area
Projection.
Map_Set, -2.99379, -75.0282, /Albers, $
XMargin=0, YMargin=0, STANDARD_PARALLEL=[-19, 21], $
LIMIT=[-3.0347, -97.5755, 12.7686, -75.0291, $
-3.0343, -52.5294, -19.1762, -75.0302]
warp = Map_Patch(peruimage, peru_lon, peru_lat)
All well and good.
Now, I wish to save this warped image, which is in a map
projection, into a GeoTiff file so I can display it in ENVI. :-)
I believe I know how to do this. I just need to project
the lat/lon values I have into UV space so I can create
a GeoTag structure. Easily done.
alberMap = MAP_PROJ_INIT('Albers Equal Area Conic', $
DATUM=8, $ ; WGS84
CENTER_LAT=-2.99379, $
CENTER_LON=-75.0282, $
STANDARD_PAR1=-19, $
STANDARD_PAR2=20)
uv = MAP_PROJ_FORWARD([peru_lon[0,0], peru_lon[0,861], $
peru_lon[1199, 861], peru_lon[1199, 0]], $
[peru_lat[0,0], peru_lat[0,861], $
peru_lat[1199, 861], peru_lat[1199, 0]], $
MAP_STRUCTURE=alberMap)
Here are the corners of my image in lat/lon and in UV coordinates,
clockwise from lower-left:
Lat/Lon in Creation Program
-99.3065 -19.4763
-98.2939 12.9533
-51.8151 12.9516
-50.8084 -19.4735
UV lat/lon in Creation Program
-2555579.0 -1894676.3
-2436402.2 1875559.7
2430893.3 1875344.8
2549426.7 -1894388.0
Looks good to me. So I then create the tie-points and scales
for my image:
s = Size(warp, /DIMENSIONS)
xscale = Abs(uv[0,0] - uv[0,2]) / s[0]
yscale = Abs(uv[1,0] - uv[1,2]) / s[1]
tp = [uv[0,1], uv[1,1]]
Here are the values (tie-point in upper-left, as usual):
Scales in Creation Program 4155.3935 4525.8355
Tie Point in Creation Program: -2436402.2 1875559.7
Then I create the geotag information and write the file:
geotag = { MODELPIXELSCALETAG: [xscale, yscale, 0], $
MODELTIEPOINTTAG: [ 0, 0, 0, tp, 0], $
GTMODELTYPEGEOKEY: 1, $
GTRASTERTYPEGEOKEY: 1, $
GEOGRAPHICTYPEGEOKEY: 4326, $
GEOGLINEARUNITSGEOKEY: 9001, $
GEOGANGULARUNITSGEOKEY: 9102, $
PROJECTEDCSTYPEGEOKEY: 32767, $
PROJECTIONGEOKEY: 32767, $
PROJCOORDTRANSGEOKEY: 11, $
PROJLINEARUNITSGEOKEY: 9001, $
PROJSTDPARALLEL1GEOKEY: -19.00000, $
PROJSTDPARALLEL2GEOKEY: 21.000000, $
PROJNATORIGINLONGGEOKEY: -75.0282, $
PROJNATORIGINLATGEOKEY: -2.99379, $
PROJFALSEEASTINGGEOKEY: 0, $
PROJFALSENORTHINGGEOKEY: 0 }
Write_TIFF, 'alber_peru.tif', Reverse(warp, 2), GEOTIFF=geotag
Next, I read the file, get the scales and tie point out of it.
They are the same as what I put in. And I use that information
to calculate the corner points of the image, like this.
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]
uv = [[xorigin, yorigin], [xorigin, yend], $
[xend, yend], [xend, yorigin]]
When I print these out:
Corner Coordinates in Display Program:
-2436402.2 -1894461.3
-2436402.2 1875559.7
2550070.1 1875559.7
2550070.1 -1894461.3
These are *not* the same as what I had before:
UV lat/lon in Creation Program
-2555579.0 -1894676.3
-2436402.2 1875559.7
2430893.3 1875344.8
2549426.7 -1894388.0
Clearly, the xscale and yscale values are wrong, since
the tie point (the 2nd value) is in the right place.
Why are they wrong? Subtle manipulation gets them
no closer. They are very wrong. The xscale value
needs to be something like 4059.46, rather than
4155.3935 for example.
The original image size is 1200x862. The warped
image size is 1199x832. I use the warped image
size in the scale equations, since that is the
image I am saving in the TIFF file and the one
that is in the map projection.
The one assumption I make that I cannot prove is
that the corners of the warped image have the same
lat/lon value as the corners of the unwarped GOES
image and the same center. Is this assumption valid?
If not, how can I determine the corner lat/lons of
the warped image?
(Usually, I have discovered my error if I have written
this much, but not this time, so I must be really stumped.)
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.")
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
