On Feb 19, 1:30 am, David Fanning <n...@dfanning.com> wrote:
> Ken Mankoff writes:
>> 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=B0S as the latitude of true scale
>> and 0=B0E 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.
>
> Oh, oh. This leads me to think you are using MAP_SET
> to set up your map projections, instead of the more
> accurate MAP_PROJ_INIT. After my mapping "epiphany" of
> a couple of weeks ago, I have given up on the MAP_SET
> projections entirely, except--possibly--for cartoon
> maps.
>
> And, of course, there is apparently a bug in the
> "more accurate" MAP_PROJ_INIT routines, in that if
> you use the UVBOX information in the map structure coming
> directly from MAP_PROJ_INIT to set up your "data
> coordinate space" for map overlays, you will still
> be "slightly off". You need to use the UVBOX information
> coming from MAP_PROJ_IMAGE for completely accurate
> results. I've had a call into ITTVIS for three weeks
> about this, but so far without results.
>
> http://www.dfanning.com/map_tips/tiffoverlay.html
>
> There could also be some confusion about whether the
> reported corner pixel coordinates are in the center
> of the pixel (likely) or on the edge of the pixel.
> If it is the center, then you are going to have to
> move the coordinates to the edge of the pixel, which
> is what IDL needs.
>
> http://www.dfanning.com/map_tips/precipmap.html
>
> Let's see what you are doing. And can your provide
> a link to an image?
>
> Cheers,
>
> David
Hi David,
Thanks for the offer to help. I went over your pages again and began
coding from scratch but am in the same place. I realize there in on
piece of info I have (71°S as the latitude of true scale) that I am
not using.
Here is the result:
http://edgcm.columbia.edu/~mankoff/tmp/unroll.png
And the code looks like this:
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=-90)
lonlat = MAP_PROJ_INVERSE( xx, yy, MAP_STRUCTURE=stereo )
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, -1*x0, -1*y0 ]
warp = MAP_PROJ_IMAGE( data, range, $
image_structure= stereo, $ ;; input
map_structure = cyl, $ ;; output
missing = -2, $
min_value = 0, $
_EXTRA=e )
erase
tv, congrid( warp, !d.x_size, !d.y_size )
pos = [0,0,1,1]
uv_box = cyl.uv_box
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
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
