Daniel Peduzzi wrote:
>
> Thanks...that is a handy program, and I've used it before in the past.
> I'm not sure that it can be used for what I want to do, though, since
> I don't want to remap the data...only display it in its *native*
> projection with coastal boundaries.
>
> In other words, if I have 1000 scanlines of DMSP data (1465 elements
> wide), and accompanying 1465x1000 lat/lon arrays, I'd like to display
> a 1465x1000 image overlaid with coastlines.
I did something like this to check if an AVHRR pixel was over land or
not.
I'll post the code here in case you are interested. I had to create an
image which consisted of 0s and 1s corresponding to sea and land. I did
this from IDL using map_set and map_continents (filling it as color=1),
then tvrd() and saving into a format handy format (in this case, idl
save format) You could try doing it without filling, just keeping the
continent outline in a file. I then have to restore this file each time
I need to check for land. This is fairly resolution dependent, but works
OK for me (AVHRR data around Antarctica).
It's a rather quick and dirty method - but *it works for me*(TM)
Ben
;=========================================================== ================
function landmaskcheck, lats, longs
; land_mask.idl is an image which has a 0 over the sea and a 1 over
; land. This was produced from a 2048x2048 window and using the
; map_set and map_continents procedures to define areas of land and sea.
; Then, using the convert_coord function we convert latitudes and
; longitudes to land mask subscripts to determine if that pixel is over
; land.
sizeimg = size(lats)
; this file contains 0s and 1s corresponding to sea/land.
; restoring this file creates an IDL variable called MASK
restore,file='~/cloudcl/data/land_mask.idl'
; open up a new window
oldwin = !D.window
window,xs = 2048,ys=2048,/pix,/free
newwin = !d.window
; setup the map reprojection used to create the land mask file initially
map_set,-90,0,0,/ster,/noborder,xmarg=0,ymarg=0,limit=[-30,- 90,-45,0,-80,90,-55,-135],/iso
; convert the input image longitudes and latitudes to device coordinates
in
; the new reprojection
sub = convert_coord(longs,lats,/data,/to_device)
; squish them to the right size
xsub = reform(sub[0,*],sizeimg[1],sizeimg[2])
ysub = reform(sub[1,*],sizeimg[1],sizeimg[2])
; this bit produces an image (same size as the input) which contains a 1
over
; land and a 0 over the ocean
flag = mask[xsub,ysub]
wdelete,newwin
wset,oldwin
return,flag
end
;=========================================================== =================
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
