| Inverse Map Projection Help [message #58811] |
Mon, 18 February 2008 17:47  |
mankoff
Messages: 131
Registered: March 2004
|
Senior Member |
|
|
Hi Group,
I've read the previous posts on inverse map projections and the
lengthy tutorial by David Fanning, but still cannot get things to line
up quite right. So I'm posting here for help...
I have a data set (BEDMAP) with this information in the header:
ncols 1371
nrows 1371
xllcorner -3426225.75
yllcorner -3426225.75
cellsize 5000
NODATA_value -9999
And this information on the website:
Polar Stereographic projection with 71°S as the latitude of true scale
and 0°E as the central meridian.
I've managed to load the data, and inverse project it approximately
such that things roughly line up. But I cannot get it accurate where
my reference for 'accurate' is the /MAP_CONTINENTS, /HIRES keywords.
My goal is to inverse project it to an equilateral lat/lon grid
(something like what comes from /CYLINDRICAL) so I can generate an
image that would work in Google Earth.
I'm pleased to post the code I've written if necessary. If anyone has
any suggestions or familiarity with this data set I would appreciate
any tips.
|
|
|
|
| Re: Inverse Map Projection Help [message #58825 is a reply to message #58811] |
Mon, 25 February 2008 06:56  |
mankoff
Messages: 131
Registered: March 2004
|
Senior Member |
|
|
On Feb 24, 2:54 pm, Paul Levine <paul.lev...@ucla.edu> wrote:
> On 2008-02-22 09:48:38 -0800, mankoff <mank...@gmail.com> said:
>
>
>
>>>> Polar Stereographic projection with 71=B0S as the latitude of true scale
>>>> and 0=B0E as the central meridian.
>
>> ;; this is the projection the data is distributed on
>> stereo = map_proj_init('Polar Stereographic', /GCTP, DATUM=8, $
>> CENTER_LONGITUDE=0, CENTER_LATITUDE=-90)
>
> You must change the CENTER_LATITUDE to -71
>
> Polar stereographic projections are free of aereal distortion at only
> one latitude, with aereal distortion increasing both north and south of
> this latitude. In the case of your data, that latitude is 71 south.
> FWIW, the northern-hemisphere sea ice data distributed by the NSIDC is
> centered at 70 north.
>
> HTH,
> Paul
It does help. Image is better aligned. But still not accurate :(.
|
|
|
|
| Re: Inverse Map Projection Help [message #58830 is a reply to message #58811] |
Sun, 24 February 2008 11:54 |
Paul Levine
Messages: 29
Registered: February 2008
|
Junior Member |
|
|
On 2008-02-22 09:48:38 -0800, mankoff <mankoff@gmail.com> said:
>>
>>> Polar Stereographic projection with 71=B0S as the latitude of true scale
>>> and 0=B0E as the central meridian.
>
>
> ;; this is the projection the data is distributed on
> stereo = map_proj_init('Polar Stereographic', /GCTP, DATUM=8, $
> CENTER_LONGITUDE=0, CENTER_LATITUDE=-90)
You must change the CENTER_LATITUDE to -71
Polar stereographic projections are free of aereal distortion at only
one latitude, with aereal distortion increasing both north and south of
this latitude. In the case of your data, that latitude is 71 south.
FWIW, the northern-hemisphere sea ice data distributed by the NSIDC is
centered at 70 north.
HTH,
Paul
|
|
|
|
|
|
| Re: Inverse Map Projection Help [message #58917 is a reply to message #58825] |
Mon, 25 February 2008 12:11 |
Paul Levine
Messages: 29
Registered: February 2008
|
Junior Member |
|
|
On 2008-02-25 06:56:50 -0800, mankoff <mankoff@gmail.com> said:
> On Feb 24, 2:54�pm, Paul Levine <paul.lev...@ucla.edu> wrote:
>> On 2008-02-22 09:48:38 -0800, mankoff <mank...@gmail.com> said:
>>
>>
>>
>>>> > Polar Stereographic projection with 71=B0S as the latitude of true s
> cale
>>>> > and 0=B0E as the central meridian.
>>
>>> ;; this is the projection the data is distributed on
>>> stereo = map_proj_init('Polar Stereographic', /GCTP, DATUM=8, $
>>> � � � � � � � � � � � �CENTER_LONGITUDE=0, CEN
> TER_LATITUDE=-90)
>>
>> You must change the CENTER_LATITUDE to -71
>>
>> Polar stereographic projections are free of aereal distortion at only
>> one latitude, with aereal distortion increasing both north and south of
>> this latitude. �In the case of your data, that latitude is 71 south. �
>
>> FWIW, the northern-hemisphere sea ice data distributed by the NSIDC is
>> centered at 70 north.
>>
>> HTH,
>> Paul
>
> It does help. Image is better aligned. But still not accurate :(.
Is the inaccuracy greater or lesser than one pixel?
|
|
|
|
| Re: Inverse Map Projection Help [message #58928 is a reply to message #58811] |
Wed, 27 February 2008 16:09 |
mankoff
Messages: 131
Registered: March 2004
|
Senior Member |
|
|
Thank you! Your second article was about shifting by half a cell
width, so I am not surprised to see this solution. I attempted that
but must have shifted in the wrong direction.
It appears like it lines up in IDL, although exporting that same image
and then importing to an overlay in Google Earth creates new alignment
problems.
Thanks again,
-k.
|
|
|
|
| Re: Inverse Map Projection Help [message #58960 is a reply to message #58811] |
Tue, 26 February 2008 20:56 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
mankoff writes:
> And the actual one I've started with is:
> http://www.antarctica.ac.uk/bas_research/data/access/bedmap/ download/surface=
> .asc.gz
Well, I just basically moved the corners over to the image
edges, which is what IDL requires, and I get what I think
is a pretty darn good fit. Here is the code I used:
pro load_asc, file, data0, data1, img
if not keyword_set(file) then begin
print, 'bathy, bedelev, groundbed, icethic, surface, water'
return
end
result = read_ascii(file+'.asc',data_start= 6)
data0 = result.field0001
bad = where( data0 eq -9999, complement= good )
data1 = data0 & data1[bad] = !values.f_nan
mm = [Min(data0[good]) , Max(data0[good])]
img = bytscl( data1, min= mm[0], max= mm[1], top= 253 ) + 1
img[ bad ] = 0
end
load_asc, 'surface', d1, d2, data
data = reverse(data,2)
s = Size(data, /Dimensions)
x0 = -2713600 -2500 ; from data set header
y0 = -2304000 -2500
x1 = s[0]*5000 + x0 + 2500
y1 = s[1]*5000 + y0 + 2500
xx = [x0,x0,x1,x1] ; the four corners
yy = [y0,y1,y1,y0]
;; this is the projection the data is distributed on
stereo = map_proj_init(106, DATUM= 8, $
CENTER_LONGITUDE= 0, CENTER_LATITUDE= -71 )
lonlat = MAP_PROJ_INVERSE( xx, yy, MAP_STRUCTURE= stereo )
longitude = reform(lonlat[0,*])
latitude = reform(lonlat[1,*])
;; output zoom
limit = [ -90, -180, max(latitude), 180 ]
;; this is the projection I would like it on
cyl = map_proj_init('Cylindrical', limit= limit)
range = [ x0, y0, x1, y1 ]
warp = MAP_PROJ_IMAGE( data, range, $
image_structure= stereo, $ ;; input
map_structure = cyl, $ ;; output
missing = -2, $
uvrange = uvrange, $
min_value = 0, $
_EXTRA= e )
erase
window, xsize=s[0], ysize=s[1]
TV, BytScl(warp)
pos = [0,0,1,1]
;; Pick one. Which one?
;uv_box = cyl.uv_box
uv_box = uvrange
Plot, uv_box[[0, 2]], uv_box[[1, 3]], Position= pos, $
/Nodata, XStyle= 5, YStyle= 5, /NoErase
MAP_CONTINENTS, Map_Structure= cyl, /HIRES
map_grid, glinest= 0, color= 255, /label, map_structure= cyl
end
--
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.")
|
|
|
|
| Re: Inverse Map Projection Help [message #58961 is a reply to message #58909] |
Tue, 26 February 2008 18:12 |
mankoff
Messages: 131
Registered: March 2004
|
Senior Member |
|
|
On Feb 26, 7:43 am, David Fanning <n...@dfanning.com> wrote:
> mankoff writes:
>> It does help. Image is better aligned. But still not accurate :(.
>
> I just got back to my office and I'm doing the usual
> up-at-3AM-thing for a week or so. Do you mean "not accurate"
> in the way using the UV-BOX from the map structure, rather
> from the MAP_PROJ_IMAGE UV-BOX, is not accurate? This wouldn't
> surprise me. Did you try using MAP_PROJ_IMAGE for creating
> the UV-BOX, as I outlined in my article?
>
> If you make the data available, I'll schedule an appointment
> for tomorrow at 4AM. :-)
>
> 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.")
Here is the data set website:
http://www.antarctica.ac.uk/bas_research/data/access/bedmap/ download/
And the actual one I've started with is:
http://www.antarctica.ac.uk/bas_research/data/access/bedmap/ download/surface.asc.gz
Code to read in this file (once un-gzipped) is:
pro load_asc, file, data0, data1, img
if not keyword_set(file) then begin
print, 'bathy, bedelev, groundbed, icethic, surface, water'
return
end
result = read_ascii(file+'.asc',data_start=6)
data0 = result.field0001
bad = where( data0 eq -9999, complement=good )
data1 = data0 & data1[bad] = !values.f_nan
mm = minmax( data0[good] )
img = bytscl( data1, min=mm[0], max=mm[1], top=253 ) + 1
img[ bad ] = 0
end
My code to attempt to 'unroll' this data is above in this thread, and
re-pasted here (slightly different than above perhaps... 2 days
later). Note that I have uv_box from both map_proj_init and
map_proj_image. I think the map_proj_image code provides slightly
better match. It appears to mach East/West perfectly (?) but there is
still a north/south error.
pro unroll_foo
end
;; load the data
load_asc,'surface', data & save, data
;restore
data = reverse(data,2)
x0 = -2713600 ; from data set header
y0 = -2304000
xx = [x0,x0,-1*x0,-1*x0] ; the four corners
yy = [y0,-1*y0,-1*y0,y0]
;; this is the projection the data is distributed on
stereo = map_proj_init('Polar Stereographic', /GCTP, DATUM=8, $
CENTER_LONGITUDE=0, CENTER_LATITUDE=-71 )
lonlat = MAP_PROJ_INVERSE( xx, yy, MAP_STRUCTURE=stereo )
longitude = reform(lonlat[0,*])
latitude = reform(lonlat[1,*])
;; output zoom
limit = [ -90, -180, max(latitude), 180 ]
;limit = [ -80, 150, -70, 180 ]
;; this is the projection I would like it on
cyl = map_proj_init('Cylindrical', limit=limit)
range = [ x0, y0, -1*x0, -1*y0 ]
warp = MAP_PROJ_IMAGE( data, range, $
image_structure= stereo, $ ;; input
map_structure = cyl, $ ;; output
missing = -2, $
uvrange = uvrange, $
min_value = 0, $
_EXTRA=e )
erase
tv, congrid( warp, !d.x_size, !d.y_size )
pos = [0,0,1,1]
;; Pick one. Which one?
uv_box = cyl.uv_box
uv_box = uvrange
Plot, uv_box[[0, 2]], uv_box[[1, 3]], Position=pos, $
/Nodata, XStyle=5, YStyle=5, /NoErase
MAP_CONTINENTS, Map_Structure=cyl, /HIRES
map_grid, glinest=0, color=255, /label, map_structure=cyl
end
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
