| Re: Convert planetary sinusoidal map image to (lat,lon,value) array [message #31626] |
Tue, 30 July 2002 18:15 |
JD Smith
Messages: 850
Registered: December 1999
|
Senior Member |
|
|
On Mon, 29 Jul 2002 11:27:08 -0700, David Fanning wrote:
> Ken Mankoff (mankoff@I.HATE.SPAM.cs.colorado.edu) writes:
>
>> I guess I'll give this a shot. Never done this but here is what I would
>> try first if I were doing it:
>
> Not too bad, not too bad. :-)
>
> It would certainly work something like this, *if* you could get an IDL
> map projection to work. I have my doubts about what a "sinusoidal equal
> area" projection is, but assuming it is the same thing as a sinusoidal
> projection in IDL, I would make a few modifications to your suggestions.
>
>> 1) Set up an identical map projection. The key word (no pun intended)
>> here is *identical*, so that every pixel of your new "sinusoidal equal
>> area" projection is the same as the BYTARR you have been given. If you
>> can't get it identical, I am not sure what to do next.
>>
>> You can test if it is identical via:
>> IDL> WINDOW, XSIZE=1440, YSIZE=720
>> IDL> MAP_SET, 0, 0, /SIN, /ISO
>> IDL> TV, sinusoidal_projection
>> IDL> MAP_GRID & MAP_CONTINENTS & MAP_HORIZON
>>
>> Does the border created by MAP_HORIZON line up exactly (to the pixel)
>> with the data your data from the TV command? Look up the keywords to
>> MAP_SET if not, and try other projections...
>
> I would add MARGIN=0 and NOBORDER=1 keywords to the MAP_SET command. And
> if I knew (or could figure out) the lat/lon coordinates of the corners
> of the image I would add the 8-element version of the LIMIT keyword as
> well. Then there is a reasonably good chance the map might match the
> image.
>
> You could set up the map coordinates in a pixmap if you didn't want to
> see something happening on the display. Just make the pixmap the same
> size as your image, etc.
>
>> 2) Redo the above code, but stop after the TV command so its just your
>> data, nothing extra. You don't actually need a window, but you need the
>> MAP_SET command to be run (with the correct 1440,720 sizes) so that IDL
>> defines the map coordinate system. You are going to use this later to
>> convert between (x,y) pixels and (lat,lon) degrees.
>>
>> 3) Set up the new array you want to put your data into. I suggest a
>> cylindrical 'projection', which can then be warped to any other
>> projection you want. So...
>> IDL> new_array = BYTARR( 1440, 720 )
>
> I would just leave it in its current projection, assuming this is the
> correct one.
>
>> 4) Step through every (x,y) pixel in your image.
>
> Since the original poster indicated that he was only interested in
> non-zero values, I would just work with those.
>
>> FOR x=0,1439 DO BEGIN
>> FOR y=0,719 DO BEGIN
>> aPixel = sinusoidal( x, y )
>> IF ( aPixel NE 0 ) THEN BEGIN ; not a valid (lat,lon) coord.
>> latlon = CONVERT_COORD( x, y, /device, /to_data ) new_array[
>> latlon[0], latlon[1] ] = aPixel
>> ENDIF
>> ENDFOR
>> ENDFOR
>
> I would do it something like this:
>
> indices = Where(sinusoidal_projection GT 0, count) IF count EQ 0 THEN
> Message, 'Whoops. Somethin gone wrong!!' x = indices MOD 1439 y =
> indices / 1439
I think you meant 1440.
> latlon = CONVERT_COORD( x, y, /DEVICE, /TO_DATA )
>
> Now, the information you want is in latlon.
Not to oversimplify, but the sinusoidal projection is given by:
x = longitude*cos(latitude)
y = latitude
You just need the inverse of this simple transform.
The y transform is trivial: knowing the range of latitude present in your
grid (lat_high and lat_low), you can write:
IDL> lat=findgen(720)/719*(lat_high-lat_low)+lat_low
OK, now you have the latitude for each row in your image. Given the known
range of longitudes present in your grid (lon_low to lon_high), you can
calculate:
IDL> lon=(findgen(1440,720) mod 1440/1439*(lon_high-lon_low)+lon_low)/ $
rebin(reform(cos(lat),1,720),1440,720)
Here I've simply linearly mapped the full range of x to the full range of
longitude (which is only valid at the equator), then divided by a suitably
inflated cos(lat). This will of course contain spurious values for the
longitude in the corners, which you can clip with:
IDL> lon[where(lon gt lon_high OR lon lt lon_low)]=!VALUES.F_NAN
And you might convince yourself it's right with:
IDL> WINDOW, XSIZE=1440, YSIZE=720
IDL> tmp=lon & tmp[where(finite(tmp) eq 0b)]=0. & tvscl,tmp
Try this with:
IDL> lat_low=-!PI/2 & lat_high=!PI/2 & lon_low=-!PI & lon_high=!PI
and you'll see a lovely map of longitude over the full globe. This also
reveals that the other suggested method (convert_coord) fails, without the
all-important POSITION keyword:
IDL> map_set,0,0,/SIN,/ISO,/NOBORDER,/NOERASE,POSITION=[0.,0.,1., 1.]
IDL> map_continents
IDL> map_horizon
IDL> map_grid,LONDEL=15,LATDEL=1
Remember that near the poles, cos=0, and the longitude is singular. What
this means in practice is that in these regions the longitude mapping is
much less certain, and subject to a greater degree of roundoff-induced
error, etc. The lat-lon bins are also severely distorted.
Good luck,
JD
|
|
|
|
|
|
| Re: Convert planetary sinusoidal map image to (lat,lon,value) array [message #31641 is a reply to message #31640] |
Mon, 29 July 2002 11:27  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Ken Mankoff (mankoff@I.HATE.SPAM.cs.colorado.edu) writes:
> I guess I'll give this a shot. Never done this but here is what I
> would try first if I were doing it:
Not too bad, not too bad. :-)
It would certainly work something like this, *if* you could
get an IDL map projection to work. I have my doubts about
what a "sinusoidal equal area" projection is, but assuming it
is the same thing as a sinusoidal projection in IDL, I would
make a few modifications to your suggestions.
> 1) Set up an identical map projection. The key word (no pun intended)
> here is *identical*, so that every pixel of your new "sinusoidal equal
> area" projection is the same as the BYTARR you have been given. If
> you can't get it identical, I am not sure what to do next.
>
> You can test if it is identical via:
> IDL> WINDOW, XSIZE=1440, YSIZE=720
> IDL> MAP_SET, 0, 0, /SIN, /ISO
> IDL> TV, sinusoidal_projection
> IDL> MAP_GRID & MAP_CONTINENTS & MAP_HORIZON
>
> Does the border created by MAP_HORIZON line up exactly (to the pixel)
> with the data your data from the TV command? Look up the keywords to
> MAP_SET if not, and try other projections...
I would add MARGIN=0 and NOBORDER=1 keywords to the MAP_SET
command. And if I knew (or could figure out) the lat/lon coordinates
of the corners of the image I would add the 8-element version of the
LIMIT keyword as well. Then there is a reasonably good chance the
map might match the image.
You could set up the map coordinates in a pixmap if you didn't want
to see something happening on the display. Just make the pixmap the
same size as your image, etc.
> 2) Redo the above code, but stop after the TV command so its just your
> data, nothing extra. You don't actually need a window, but you need
> the MAP_SET command to be run (with the correct 1440,720 sizes) so
> that IDL defines the map coordinate system. You are going to use this
> later to convert between (x,y) pixels and (lat,lon) degrees.
>
> 3) Set up the new array you want to put your data into. I suggest a
> cylindrical 'projection', which can then be warped to any other
> projection you want. So...
> IDL> new_array = BYTARR( 1440, 720 )
I would just leave it in its current projection, assuming
this is the correct one.
> 4) Step through every (x,y) pixel in your image.
Since the original poster indicated that he was only
interested in non-zero values, I would just work with
those.
> FOR x=0,1439 DO BEGIN
> FOR y=0,719 DO BEGIN
> aPixel = sinusoidal( x, y )
> IF ( aPixel NE 0 ) THEN BEGIN ; not a valid (lat,lon) coord.
> latlon = CONVERT_COORD( x, y, /device, /to_data )
> new_array[ latlon[0], latlon[1] ] = aPixel
> ENDIF
> ENDFOR
> ENDFOR
I would do it something like this:
indices = Where(sinusoidal_projection GT 0, count)
IF count EQ 0 THEN Message, 'Whoops. Somethin gone wrong!!'
x = indices MOD 1439
y = indices / 1439
latlon = CONVERT_COORD( x, y, /DEVICE, /TO_DATA )
Now, the information you want is in latlon.
Cheers,
David
--
David W. Fanning, Ph.D.
Fanning Software Consulting, Inc.
Phone: 970-221-0438, E-mail: david@dfanning.com
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Toll-Free IDL Book Orders: 1-888-461-0155
|
|
|
|
| Re: Convert planetary sinusoidal map image to (lat,lon,value) array [message #31642 is a reply to message #31641] |
Mon, 29 July 2002 10:58 |
Ken Mankoff
Messages: 158
Registered: February 2000
|
Senior Member |
|
|
Hi Suniti,
I guess I'll give this a shot. Never done this but here is what I
would try first if I were doing it:
1) Set up an identical map projection. The key word (no pun intended)
here is *identical*, so that every pixel of your new "sinusoidal equal
area" projection is the same as the BYTARR you have been given. If
you can't get it identical, I am not sure what to do next.
You can test if it is identical via:
IDL> WINDOW, XSIZE=1440, YSIZE=720
IDL> MAP_SET, 0, 0, /SIN, /ISO
IDL> TV, sinusoidal_projection
IDL> MAP_GRID & MAP_CONTINENTS & MAP_HORIZON
Does the border created by MAP_HORIZON line up exactly (to the pixel)
with the data your data from the TV command? Look up the keywords to
MAP_SET if not, and try other projections...
2) Redo the above code, but stop after the TV command so its just your
data, nothing extra. You don't actually need a window, but you need
the MAP_SET command to be run (with the correct 1440,720 sizes) so
that IDL defines the map coordinate system. You are going to use this
later to convert between (x,y) pixels and (lat,lon) degrees.
3) Set up the new array you want to put your data into. I suggest a
cylindrical 'projection', which can then be warped to any other
projection you want. So...
IDL> new_array = BYTARR( 1440, 720 )
4) Step through every (x,y) pixel in your image. Actually, don't use
the image itself, just use your array. For each (x,y), figure out what
(lat,lon) it is at. Stick it into the correct bin in NEW_ARRAY.
FOR x=0,1439 DO BEGIN
FOR y=0,719 DO BEGIN
aPixel = sinusoidal( x, y )
IF ( aPixel NE 0 ) THEN BEGIN ; not a valid (lat,lon) coord.
latlon = CONVERT_COORD( x, y, /device, /to_data )
new_array[ latlon[0], latlon[1] ] = aPixel
ENDIF
ENDFOR
ENDFOR
I am not sure if my syntax is correct, but this should give you the
general idea
Hope this helps,
Let us know if it works,
-Ken.
On 29 Jul 2002, Suniti Karunatillake wrote:
> I have a sinusoidal equal area map projection image in the form of a
> byte array. The image is rectangular, 1440 pixels (length) X 720
> pixels (height). I am able to read the file and store it as an array
> of dimensions 1440(columns)X720(rows).
>
> I need to convert the planetary portion of the image (all other points
> have value 0 in the image array) into a table of form (latitude,
> longitude, value). I am also uncertain which type of array I should
> construct from the planetary portion of the image in order to apply
> the conversion.
>
> Sincerely,
> Suniti
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
