IDL projections (MAP_PROJ_IMAGE) and ENVI projections, Select spatial substets of images

[sh]

Hi together,

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??

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)

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??

- I use congrid to resize the image to a "plotable" size. Maybe this
causes the shifting between the coastlines?

- 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?

- The images are correct but the coastlines are in a wrong
projection??

Maybe some of you have already expierence with such tasks! Any help
would be appreciated!


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)

window, /free, xsize=10983./10., ysize=7797./10.
disp=congrid(warped4, 10983./10., 7797./10.)
tv,reverse(disp,2)
plot, map4.uv_box[[0,2]],map4.uv_box[[1,3]],/NoData, XSTYLE=1,
YSTYLE=1,POSITION=[0,0,1,1], /noerase
map_continents, map_structure=map4
map_grid, map_structure=map4, londel=5, latdel=5


; to display ENVI result
infile = 'D:\temp\test.img'
read_envi_file, infile, img, xs, ys, type, offset, mapinfo,
STATUS=status

window, /free, xsize=10983./10., ysize=7797./10.
disp=congrid(img, 10983./10., 7797./10.)
tv,reverse(disp,2)
map_set, /mercator, limit=geographical_extend, /noerase, /
hires,xmargin=0, ymargin=0
map_continents, /coasts, color=255, /hires
map_grid, londel=5, latdel=5



cheers,
Sebastian

[David Fanning]

Post by sh
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.

Post by sh
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.)

Post by sh
- 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.

Post by sh
- 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.

Post by sh
- 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.

Post by sh
; 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

[sh]

Post by David Fanning

Post by sh
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.

Post by sh
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.)

Post by sh
- 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
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.

Post by sh
- 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.

Post by sh
- 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.

Post by sh
; 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
   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.")

Thanks a lot David for your detailed answer!

I already saw your example of how to warpe images, but this only shows
how to project already projected images. And my image is "raw"
lat,lon.

Anyway, at least I found a (very old) program which does and can
exactly what I want. Maybe someone could take a look, if the results
are correct!? Or where is room for improvement?!

http://cimss.ssec.wisc.edu/~gumley/imagemap.html
http://cimss.ssec.wisc.edu/~gumley/idl/imagemap.pro

thx!

[David Fanning]

Post by sh
I already saw your example of how to warpe images, but this only shows
how to project already projected images. And my image is "raw"
lat,lon.

I'm not sure what "raw" means. But these are
probably just geographic coordinates. This is
sometimes called a plate carree map projection,
or--in GCTP parlance--an equirectagular projection.

Post by sh
Anyway, at least I found a (very old) program which does and can
exactly what I want. Maybe someone could take a look, if the results
are correct!? Or where is room for improvement?!

"Correct" in what sense? It is certainly not the
kind of projected data you can do science with, as
Liam points out in his notes. It is a "display"
solution. In this sense, if it "looks good" it is
probably correct enough.

Cheers,

David

[Dave Poreh]

Post by David Fanning

Post by sh
I already saw your example of how to warpe images, but this only shows
how to project already projected images. And my image is "raw"
lat,lon.

I'm not sure what "raw" means. But these are
probably just geographic coordinates. This is
sometimes called a plate carree map projection,
or--in GCTP parlance--an equirectagular projection.

Post by sh
Anyway, at least I found a (very old) program which does and can
exactly what I want. Maybe someone could take a look, if the results
are correct!? Or where is room for improvement?!

"Correct" in what sense? It is certainly not the
kind of projected data you can do science with, as
Liam points out in his notes. It is a "display"
solution. In this sense, if it "looks good" it is
probably correct enough.
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.")

As Liam point out, this routine is not working for high latitudes
(weird) too!!!

[sh]

Post by Dave Poreh

Post by David Fanning

Post by sh
I already saw your example of how to warpe images, but this only shows
how to project already projected images. And my image is "raw"
lat,lon.

I'm not sure what "raw" means. But these are
probably just geographic coordinates. This is
sometimes called a plate carree map projection,
or--in GCTP parlance--an equirectagular projection.

Post by sh
Anyway, at least I found a (very old) program which does and can
exactly what I want. Maybe someone could take a look, if the results
are correct!? Or where is room for improvement?!

"Correct" in what sense? It is certainly not the
kind of projected data you can do science with, as
Liam points out in his notes. It is a "display"
solution. In this sense, if it "looks good" it is
probably correct enough.
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.")

As Liam point out, this routine is not working for high latitudes
(weird) too!!!

Finally it works like your example on your webpage David. (http://
www.dfanning.com/map_tips/warpimage.html)

The image was flipped horizontally and I just changed this at TV not
before I'm reprojecting the image!!
Stupid mistake!

[David Fanning]

Post by sh
The image was flipped horizontally and I just changed this at TV not
before I'm reprojecting the image!!

Ah, yes. Flipping images (I always use REVERSE rather than
!ORDER or absolute chaos ensues) is just something I've
learned to do always. I don't even think about it any more,
and often forget I have done it. Mapping *always* means the
(0,0) location is in the upper-left corner. This is just
something you have to memorize in IDL, like COLUMN-row.

Cheers,

David