Re: [ferret_users] Finding the area a certain distance/number of grid cells from a boundary

[Rob Briggs <rdbriggs@xxxxxx>]

Russell

Perfect and elegant. Thank you kindly. 

Rob 

P.S. This was my ham-fisted un-ferret-like solution I came up with prior
to your nice solution. I have a lot to learn before I can think like a
ferret. I wish I'd know about fill_xy previously. 

let sr = 2 ! search radius
let ocean = 0 
let start = 1
let finish = 141

! now find ice shelf/ocean boundary
REPEAT/RANGE=`start+sr`:`finish-sr`:1/name=r (\
 REPEAT/RANGE=`start+sr`:`finish-sr`:1/name=s (\
  let iob = IF aHShf[I=`r`,J=`s`] EQ 1 AND\
  (((aHSHF[I=`r+sr`,J=`s-sr`]) eq `ocean` ) OR\
   ((aHSHF[I=`r+sr`,J=`s`   ]) eq `ocean` ) OR\
   ((aHSHF[I=`r+sr`,J=`s+sr`]) eq `ocean` ) OR\
   ((aHSHF[I=`r`   ,J=`s-sr`]) eq `ocean` ) OR\
   ((aHSHF[I=`r`   ,J=`s+sr`]) eq `ocean` ) OR\
   ((aHSHF[I=`r-sr`,J=`s-sr`]) eq `ocean` ) OR\
   ((aHSHF[I=`r-sr`,J=`s`   ]) eq `ocean` ) OR\
   ((aHSHF[I=`r-sr`,J=`s+sr`]) eq `ocean` ))\
 THEN 1 ELSE 0;\
  LIST `r`,`s`,`iob`;\
  LIST/file=iob`sr`.dat/append/nohead iob;\
 )\
)




On Mon, 2011-02-21 at 13:19 +1100, Russell Fiedler wrote:
> Rob,
> 
> I think constructing a mask using the FILL_XY() function is what you want.
> 
> You have the array thick
> 
> let ocean_mask = if ( thick eq 0 ) then 1
> let shelf_mask = if ( thick gt 0 ) then 1
> let extended_ocean_mask = fill_xy(ocean_mask,shelf_mask,2)       ! All shelf points within 2 of the ocean are now also filled with 1
> 
> let shelf_near_ocean =  if ( thick gt 0 ) then extended_ocean_mask  !  Only want shelf points not ocean
> 
> shade thick
> shade/ov/pat=tiny_grid/pal=black shelf_near_ocean   ! Check
> 
> Multiplying by the area of the cells (1600 in your case) and summing should do the trick
> 
> let area = 1600 * shelf_near_ocean[i=@sum,j=@sum]
> 
> This could also be done using the @SHF operator but it's rather messy.
> 
> Russ
> 
> 
> 
> On Sunday 20 February 2011 04:12, Rob Briggs wrote:
> > Hello, 
> > 
> > I have a data set that is ice shelf thickness of Antarctica. I want to
> > find the total area of shelf that is bounded between the shelf edge and
> > say 100km inland. 
> > 
> > The elements of the dataset are either: 
> > a) some value for iceshelf thickness 
> > b) 0 for ocean 
> > c) NaN for not ice shelf and not ocean. 
> > 
> > See attached pdf. Thus there is a way to identify the ice-shelf edge,
> > i.e. boundary between ocean and ice shelf. 
> > 
> > An equivalent example, using a tweaked etopo (except my dataset is on a
> > polar stereo grid) would be, find the area of land that <500m elevation
> > and is within 100km of the ocean/land interface. 
> > 
> > yes? use etopo60
> > yes? let r1 = if rose gt 0 then rose else 0
> > yes? let r2 = if r1 lt 500 then r1
> > yes? shade r2
> > 
> > As my data is on a polar stereo grid, it doesn't even need to be a 100km
> > segment; each grid cell is know dimension e.g. 40km x 40km so if I could
> > identify the first 2 or 3 grid cells from the shelf edge grid cell that
> > would suffice.
> > 
> > Regards,
> > 
> > Rob 
> > 
> > -----------------------------------------------
> > Ph.D candidate Glacial Systems Dynamics
> > 
> > C4043C
> > Department of Physics and Physical Oceanography
> > Memorial University of Newfoundland, St. John's
> > NL, A1B 3X7, Canada
> > 
> > Phone   : 001 (709) 864-2407
> > Email   : rdbriggs@xxxxxx
> > Website : http://www.physics.mun.ca/~rdbriggs/
> >