comp.lang.idl-pvwave archive: archive » Inverse Map Projection Help

Home » Public Forums » archive » Inverse Map Projection Help

Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend

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.")

Report message to a moderator

[Message index]

		Inverse Map Projection Help By: mankoff on Mon, 18 February 2008 17:47
		Re: Inverse Map Projection Help By: mankoff on Mon, 25 February 2008 06:56
		Re: Inverse Map Projection Help By: David Fanning on Tue, 26 February 2008 04:43
		Re: Inverse Map Projection Help By: mankoff on Tue, 26 February 2008 18:12
		Re: Inverse Map Projection Help By: Paul Levine on Mon, 25 February 2008 12:11
		Re: Inverse Map Projection Help By: Paul Levine on Sun, 24 February 2008 11:54
		Re: Inverse Map Projection Help By: mankoff on Wed, 27 February 2008 16:09
		Re: Inverse Map Projection Help By: David Fanning on Tue, 26 February 2008 20:56

Previous Topic:	How to get file UNIT from file NAME?
Next Topic:	Re: How to get numbers into passed structure elements (pass-by-value/reference problem).

-=] Back to Top [=-

[ Syndicate this forum (XML) ] [

] [

]

Current Time: Sun Oct 12 01:58:01 PDT 2025

Total time taken to generate the page: 0.01112 seconds