[Thread Prev][Thread Next][Index]

Re: [ferret_users]plotting transect along the deepest depth (i.e thalweg)



Hi Murat,
          My lack of awareness about "thalweg" caused the questions in my
previous mail. Thanks to Google and Wikipedia. For those who are interested
to know about "thalweg" here is the link : http://en.wikipedia.org/wiki/Thalweg

  Using Levitus dataset, it is difficult to create an example solution for 
your question. Hence, I have flipped the "valley"s in your question to "ridges".
The procedure remains the same (especially for finding deepest depth).
But, when you deal with your actual case ("thalweg"), you have to change
 
             let z_crit       = topo_xymask[i=@MIN]
to 
             let z_crit       = topo_xymask[i=@MAX]

See the detailed demo below and try to adapt this demo for your case....
Also note that you don't have to use etopo20 data set, if your variable
is on a high resolution Z-axis. Please mail back to the list, if you have
any questions/problems.

Regards,

Jaison,

!-------------------------DEMO starts here------------------------------
\ cancel mode verify
!
! Description : A demo to plot the temperature values at the deepest depth
!                    along a ridge (90E ridge). It can be easily adapted to
!                    find the same along a thalweg. See the comments below.
!
! Created By : JK
! Creted  On : 19/June/2007
!
!-------------------------------------------------------------------------- 

  use levitus_climatology
  use etopo20
 
  let field = temp[d=1]
  let topo  = rose[d=2,gx=field,gy=field]*-1  ! flip the sign, and
                                              !  assign field's grid
  shade/lev=(0,6000,250)/x=30:120/y=-50:30 topo
  go land

  ! locate 90E ridge : See the X-axis values precisely

  go box 85.5,91.5,(-40.5),15.5 
 
  ! mask other regions 
  let mask_x = IF x[gx=topo] GE 86.5 AND x[gx=topo] LE 91.5 THEN 1
  let mask_y = IF Y[gy=topo] GE (-40.5) AND Y[gy=topo] LE 15.5 THEN 1
   
  let topo_xymask = topo * mask_x * mask_y

  pause ; shade /lev=(0,6000,250)/x=30:120/y=-50:30 topo_xymask
  go land

  ! locate grid points with shallowest depth (topo_xymask[i=@MIN]) in 
  !    the xy masked data to mark the grid points on a ridge.  
  !    If you are dealing with a thalweg, find the deepest depth
  !    (topo_xymask[i=@MAX]).
  ! 
  ! Then create a final mask (mask_z).
  
  let z_crit       = topo_xymask[i=@MIN]  
  let zero_at_crit = topo_xymask - z_crit
  let mask_z       = zero_at_crit[i=@WEQ:0]
  
  pause ; shade/ov /lev=(0,6000,250)/x=30:120/y=-50:30 mask_z

  ! apply the topo mask

  let field_masked = field + mask_z * 0 

  ! Find the deepest depth along ridge/thalweg. Then use it to find  
  !    the value of given field at deepest depth.
  ! 
  ! remember that the following lines will not change even if you
  !    intend to work on so called "thalweg"

  let kvals         = k[g=field_masked] + field_masked * 0
  let kzero         = kvals - kvals[k=1:`field_masked,return=kend`@MAX]
  let field_deep    = field_masked * kzero[k=@WEQ:0]
  let field_deep_1D = field_deep[k=@SUM,i=@SUM]

  pause ; shade/x=30:120/y=-50:30 field_deep[k=@SUM]

  pause ; plot/y=-50:30 field_deep_1D

!-------------------------DEMO ends here------------------------------




On Tue, 19 Jun 2007, Jaison Kurian wrote:

> Hi Murat,
>           Two questions :
>  
>    1. What exactly is your problem? (Or what was the error message or
>          difficulty that you have faced?)  
> 
>    2. While finding kzero, what the "I=1:`temperature,RETURN=iend`@max"
>          stands for? If you intend to do the a similar operation on X 
>          axis, you have to repeat all the steps for X-axis too. 
> 
> The example given in the webpage you have mentioned is working fine for
> me. The procedure is explained step by step below :
> 
>   --->  let kvals = k[g=temp] + temp * 0
> 
>    This step will create a new variable "kvals", with following features 
> 
>        - XYZT grid will be EXACTLY same as that of "temp"
>        - at each and every grid point, the value of kvals will be the
>             level number (K-index) corresponding to the grid point
>        - Missing values will be EXACTLY as in "temp"
>    
>    Example :   yes? use levitus_climatology
>                yes? let kvals = k[g=temp] + temp * 0
>                yes? list/x=80/y=10 temp, kvals
> 
>   ---> let kzero = kvals - kvals[k=1:`temp,return=kend`@MAX]
> 
>    This step will create a new variable kzero with following features
>      
>       - XYZT grid same as that of temp or kvals
>       - will have zero value along Z-axis, ONLY at the deepest depth
> 
>    Example : After the lines in above above example, issue :
> 
>               yes? kzero = kvals - kvals[k=1:`temp,return=kend`@MAX]
>               yes? list/x=80/y=10 temp, kvals, kzero
> 
>   ---> let integrand = temp * kzero[k=@WEQ:0]
> 
>     kzero[k=@WEQ:0] will return 1 if it finds a zero value exactly on
>         a point along Z-axis, else it will return the "weights" to get
>         zero, with missing values in all other Z-points. See User
>         Manual, Ch3 Sec2.4.27.  @WEQ for details.
> 
>     Example : list/x=80/y=10 temp, kvals, kzero, kzero[k=@WEQ:0]
> 
>     When you multiply temp with kzero[k=@WEQ:0], the resulting variable
>     will have valid data ONLY on those points where kzero[k=@WEQ:0] 
>     have "non-misssing" values. Now it can be directly SUMMED to get
>     the desired variable 
> 
>     Example :  yes? let integrand = temp * kzero[k=@WEQ:0]
>                yes? plot integrand[k=@SUM,x=150]
> 
> Hope this helps,
> 
> Regards,
> 
> Jaison
>           
> 
> On Tue, 19 Jun 2007, Murat Gunduz wrote:
> 
> > 
> > Dear Ferret Users,
> > 
> > I would like to make a transect plot of temperature (from south to north)
> > following the deepest depth from my 3D netcdf file.
> > 
> > I have searched the e-mail list, but I could not find any clue.
> > I tried to modify the one of the script in FAQ, but without success
> > http://ferret.wrc.noaa.gov/Ferret/FAQ/analysis/deepest_depth.html
> > 
> > 
> > my simple script is below,
> > 
> > use temperature_3z.nc
> > LET kvals = k[G=temperature] + 0*temperature
> > let kzero = kvals - 
> > kvals[K=1:`temperature,RETURN=kend`@MAX,I=1:`temperature,RETURN=iend`@max]
> > LET integrand = TEMPERATURE* kzero[k=@WEQ:0]
> > 
> > 
> > 
> > Could you please help, how can I generate 2D transect along the thalweg?
> > Thank you very much in advance.
> > 
> > Murat Gunduz
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> 
> 

-- 
This message has been scanned for viruses and
dangerous content by MailScanner, and is
believed to be clean.



[Thread Prev][Thread Next][Index]

Contact Us
Dept of Commerce / NOAA / OAR / PMEL / TMAP

Privacy Policy | Disclaimer | Accessibility Statement