[Thread Prev][Thread Next][Index]
Re: Draw coastal line in different sea depth levels.
Here is a way to draw the boundary of a mask using straight line
segments for the top, bottom, left, and right edges of grid cells
lying along the "land-sea" border.
Example: The Gulf of Alaska
use etopo20
region/x=180:270/y=45:60
shade rose
let mask=if(rose[d=1] gt 0)then 1 else 0
go bounder
where bounder.jnl is the following script
! bounder.jnl : drawing the "land-sea" edge in a mask
def grid/like=mask gm
let xl=x[g=gm]-xbox[g=gm]/2 ; let xr=x[g=gm]+xbox[g=gm]/2
let yb=y[g=gm]-ybox[g=gm]/2 ; let yt=y[g=gm]+ybox[g=gm]/2
! the left edge
let ylb=if(mask[x=@ddb] ne 0)then yb else 999
let ylt=if(mask[x=@ddb] ne 0)then yt else 999
! the right edge
let yrb=if(mask[x=@ddf] ne 0)then yb else 999
let yrt=if(mask[x=@ddf] ne 0)then yt else 999
! the bottom edge
let xbl=if(mask[y=@ddb] ne 0)then xl else 999
let xbr=if(mask[y=@ddb] ne 0)then xr else 999
! the top edge
let xtl=if(mask[y=@ddf] ne 0)then xl else 999
let xtr=if(mask[y=@ddf] ne 0)then xr else 999
! list all edges to the file bounds.d ...
list/nohead/form=(4f8.2,' 999 999')/file=bounds.d/clob xl,ylb,xl,ylt
list/nohead/form=(4f8.2,' 999 999')/file=bounds.d/app xr,yrb,xr,yrt
list/nohead/form=(4f8.2,' 999 999')/file=bounds.d/app xbl,yb,xbr,yb
list/nohead/form=(4f8.2,' 999 999')/file=bounds.d/app xtl,yt,xtr,yt
! ... then use grep to exclude the interior cells ...
sp rm bounds2.d ! this may require the user to type y
sp grep -v 999.00 bounds.d > bounds2.d
file/form=free/var=xx,yy/col=6 bounds2.d
set var/bad=999 xx ; set var/bad=999 yy
can region
plot/o/nolab/vs/line=1 xx,yy
------
Note that if a very large region were set, or the mask had a very fine
resolution, the number of points in bounds2.d might exceed the 20000
or so default limit and the read statement would need to be replaced
with a gridded read.
Mick Spillane
|--****--****-*---*---***--***--|____spillane@pmel.noaa.gov____|
|-*__---*-----*--*-*--*--*-*--*-|_SCIENCE APPLICATIONS SUPPORT_|
|--***--*-----*-*---*-***--***--|____EPIC/Ferret/PlotPlus______|
|-----*-*-----*-*****-*----*----|__Room 2070 Bldg#3 NOAA/PMEL__|
|-****---****-*-*---*-*----*----|____Phone_:_(206)526-6780_____|
[Thread Prev][Thread Next][Index]
Dept of Commerce /
NOAA /
OAR /
PMEL /
TMAP
Contact Us | Privacy Policy | Disclaimer | Accessibility Statement