Re: IDL projections (MAP_PROJ_IMAGE) and ENVI projections, Select spatial substets of images [message #71209] |
Sun, 06 June 2010 10:32 |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Sebastian writes:
> I have a problem using map projections since I didn't get the same
> results for IDL and ENVI.
>
> My task is to reproject an image in lat/lon's (WGS 84, Geographic lat
> lon) to mercator for Australia and show/save (*.png) only spatial
> defined subsets of this reprojected image. (And also define "bigger"
> spatial subset and pin the image on the right position)
>
> I tested MAP_IMAGE and MAP_PROJ_IMAGE in IDL, but in both cases a had
> a problem with the output dimension. MAP_IMAGE seems to act in
> accordance with the predefined windows size and with MAP_PROJ_IMAGE it
> is possible to set it by yourself. On the contrary, ENVI makes a
> suggestion concerning the output pixelsize, but I didn't get it how
> ENVI calculates this output size.
> So I guess the sampling rate and/or pixel size is responsible for the
> suggetion of ENVI?? Or even the distortion introduced by the map
> projection??
To transform one map projection into another you have
to put all the pixel information into XY coordinates
(sometimes called UV coordinates in IDL). In other
words, you have to do this in projected coordinate
space, not geographic coordinate (lat/lon) space.
The pixel resolution you are refering to is the pixel
spacing in this projected coordinate space. Or, another
way of saying this, the projected grid spacing. The
output resolution of the image is how many grid units
you want to have in your final image.
To do this gridding, you create a projected grid and for
every grid unit, you find the latitude and longitude
of that point (Map_Proj_Inverse). Then, you calculate
a new value by interpolating that point from your original
image values. This is basically what Map_Proj_Image does
for you.
> Anyway, for my further comparison between ENVI and IDL I used the
> output dimension suggested by ENVI to reproject the image. And to
> compare the results I made to plots with the coast lines, which in
> BOTH (!) cases didn't match (the result of ENVI was a little bit
> better somehow)
Humm. Don't know about this. I've never had trouble putting
map boundaries on images, if I have set up the map coordinate
space in the correct way. Mostly this means setting up a
projected grid range (rectangular array) so that I can "plot"
the map boundaries on it. I use the MapCoord object in my
Catalyst Library to do this. It always does a good job for
the map projections I use. (I use the Map_Outline object to
draw the map boundaries and the Map_Grid object to draw map
grids.)
> Now I have several questions/comments:
>
> - I have seen that there are 2 libraries within IDL (IDL, GCTP), so I
> tested both of them. The only difference I realized was that you can't
> set an ellipsoid for the IDL (map_set) library. Which library uses
> ENVI??
Neither. :-)
At the time ENVI was written IDL only had the Map_Set map
projections. There were (and still are!) completely inadequate
for the kind of precision ENVI wanted, so the ENVI folks wrote
their own map projection software, which I presume they still
use.
Later, IDL added the GCTP map projections and these are
much better (and the only ones you should be using if you
want professional map projections), but there have always
been problems with them (several of which are *finally*
fixed in the upcoming IDL 8) and they are also becoming
a bit long of tooth. Better open source map projections
exist (proj4 routines) and, in my mind, should be incorporated
into IDL if the folks at ITTVIS want to be current with
what's going on in this field.
> - I use congrid to resize the image to a "plotable" size. Maybe this
> causes the shifting between the coastlines?
Well, I guess it depends on how you are using it. :-)
I always use TVImage to display my images, and it uses
Congrid, of course. As I say, I've never had problems
aligning boundaries on images.
> - How can I select a spatial subset from the image and plot it into a
> "bigger" spatial subset? e.g. to show only the east coast of australia
> but with new zealand (where no image data is available) I think the
> "problem" here is to find the right position?
This is really just a gridding problem. One of the
weaknesses of IDL is that it doesn't really allow
the kind of map image gridding you need for working
with images in map projections. (I say this with
some trepidation because I am convinced that IDL
probably *does* provide this kind of support, but
in five or six years of trying to use IDL to do it,
I have come up completely empty.) At NSIDC, this
kind of gridding is done with our mapx utilities,
which are only available on UNIX platforms.
In any case, it is probably not too difficult to
produce the kind of output you want, if you set your
"data" coordinate space up correctly. Again, I would
rely on my MapCoord object to do this.
> some lines of my code:
>
> ; to display IDL result
>
> geographical_extend= [-39.5,112.5,-10.5,154.0]
> range =
> [geographical_extend[1],geographical_extend[0],geographical_ extend[3],geographical_extend[2]]
>
> ; c is the image with the size 9960, 6960 and pixelspacing 0.00417°
> map4 = map_proj_init(105, ellipsoid=8, limit=geographical_extend)
> warped4 = map_proj_image(c, range, dimensions=[10983,7797],
> map_structure=map4, uvrange=uvOut4, xindex=xindex4, yindex=yindex4)
To do this correctly, the "range" should not be in
geographical coordinates, but in projected XY coordinates.
Earlier documentation (i.e., prior to IDL 7.1) of
Map_Proj_Image was misleading on this point.
I gave a talk on IDL map projections at the last IDL
Users Group meeting which you may find helpful. You can
find my Powerpoint presentation here:
http://www.dfanning.com/powerpoint/map_projections_idl.pdf
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.")
|
|
|