[ncl-talk] How to extract data just over the circle defined by "nggcog" function

Mary Haley haley at ucar.edu
Fri Jan 19 14:35:36 MST 2018


Hi Nafiseh,

Thanks for your patience.  I was unable to respond to this question until I
had a chance to look at the documentation.

You are correct, this function is confusing.  I admit I wasn't familiar
with it until I started looking at the internal code and tried to figure
out why the restriction is placed on it.

As far as I can tell, there's no reason the two sets of arguments should
have the weird set of restrictions of dimensionality placed on them.

The first two input arrays should be allowed to be any dimensionality of
lat/lon points  that you want to loop through and check if it's on the
given lat/lon great circle, specified by the second two arguments.

I've created a trouble-ticket for this, NCL-2717, and marked it as a high
priority for NCL V6.5.0.

Meanwhile, the only work-around I can think of is to loop through each of
your points and pass them individually to gc_onarc.  This means you will
need to create the output array before you call the function, and then
subscript it as you go.  Something like this (UNTESTED, because I'm not
sure of the exact dimensionality of your arrays):

npts = dimsizes(lat_2d_outer_core_to_1d)
onarc_outer_core = new(npts,logical)

do n=0,npts-1
  onarc_outer_core(n) = gc_onarc(lat_2d_outer_core_to_
1d(n),lon_2d_outer_core_to_1d(n),clat_outer_core_2d,clon_outer_core_2d)
end do

This is going to be slower than calling the function once, but hopefully it
still works.

--Mary




On Wed, Jan 17, 2018 at 11:41 PM, nafiseh pegahfar via ncl-talk <
ncl-talk at ucar.edu> wrote:

> To whom it may concern,
> I decided to use wrf outputs and a tropical cyclone center location (from
> JTWC) to define some circular nests around TC center with fixed radiuses.
> To this aim, I used nggcog and gc_inout functions to define the circular
> nest and determine the points that are in the defined circle. Now I need to
> now which points are located on the circle. So I used gc_onarc function,
> but defining the inputs for this function confused me. the two rightmost
> inputs for gc_onarc are nggcog out puts with (/72,2/)dimension and
> dimension of the last two left gc_onarc inputs are 3600. I do highly
> appreciate if help me. The sample of code are attached as below:
>
>
>
> ;oooooooooooooooooooooooooooooooooooooooooooooooooooooooooo***************
> ;
> ;
> load "/usr/share/ncarg/nclscripts/csm/gsn_code.ncl"
> load "/usr/share/ncarg/nclscripts/csm/gsn_csm.ncl"
> load "/usr/share/ncarg/nclscripts/csm/contributed.ncl"
> load "/usr/share/ncarg/nclscripts/csm/shea_util.ncl"
> load "/usr/share/ncarg/nclscripts/wrf/WRFUserARW.ncl"
> load "/usr/share/ncarg/nclscripts/csm/skewt_func.ncl"
>
>  begin
> ;oooooooooooooooooooooooooooooooooooooooooooooooooooooooooo**************
> ; add data from  file
> ;ooooooooooooooooooooooooooooo*********
>   DATADir = "/home/pegahfar/inio/haiyan/data/wrfout_18_6_haiyan_MPI_8/"
> FILES = systemfunc (" ls -1 " + DATADir +"wrfout*")
> print ("FILES="+FILES)
>
> lat_center = (/(/5.8, 5.8, 6.1, 6.1,6.1, 6.2, 6.3, 6.5,6.5, 6.5, 6.9,
> 7.1,7.3, 7.6, 7.9, 8.2,8.7, 9.3, 10.2, 10.6,11.0, 11.4, 11.9, 12.2,12.3,
> 13.5, 14.4, 15.4,16.5, 17.9, 19.4, 20.4,21.5, 22.4, 22.4 , 22.4/)/)
> ;;;print (lat_center)
>
> lon_center =  (/(/157.2,157.2,155.5,153.3,152.2,150.4,148.8,147.2,145.9,
> 144.6,142.9,141.3,139.7,138.0,136.2,134.4,132.8,131.1,129.1,
> 126.9,124.8,122.5,120.5,118.0,116.6,114.8,113.1,111.4,110.3,
> 109.0,108.0,107.5,107.1,107.7,107.7,107.7/)/)
>
> do ifil =1,35,1
>
> a= addfile(FILES(ifil),"r") ; Open the next file
> times := wrf_user_getvar(a,"times",-1)  ; get all times in the file
> ntimes = dimsizes(times)         ; number of times in the file
>
>     lat_wrf := a->XLAT(0,:,:)
>     lon_wrf := a->XLONG(0,:,:)
>
> radius_outer_core=(7*RMW_km(ifil))/110
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;;;;;;;;;;;;;;;;;;;
> ;;;;;;;;;   /\ /\ /\ /\ start of define circle sized structure /\ /\ /\ /\
> /\;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;;;;;;;;;;;;;;;;;;;
>
>   npoint_outer_core=144
>
>   clat_outer_core = new(npoint_outer_core,float)   ; arrays to hold great
> circle
>   clon_outer_core = new(npoint_outer_core,float)
>
>     lat_location = lat_center(ifil)
>     lon_location = lon_center(ifil)
>
>     nggcog(lat_location,lon_location,radius_outer_core,
> clat_outer_core,clon_outer_core)   ; Calculate great circle
>
>     min_lat_outer_core = min(clat_outer_core)
>     min_lon_outer_core = min(clon_outer_core)
>     max_lat_outer_core = max(clat_outer_core)
>     max_lon_outer_core = max(clon_outer_core)
>
> extrem_outer_core_lat=(/min_lat_outer_core , max_lat_outer_core /)
> extrem_outer_core_lon=(/min_lon_outer_core , max_lon_outer_core /)
>
> extrem_outer_core_inds=getind_latlon2d (lat_wrf,lon_wrf,
> extrem_outer_core_lat, extrem_outer_core_lon)
>
> ;---Subset the desired rectagle of data
>
> loc_outer_core_circle:= getind_latlon2d (lat_wrf,lon_wrf, clat_outer_core,
> clon_outer_core)
>
> lat_ind_over_c=loc_outer_core_circle(:,0)
> lon_ind_over_c=loc_outer_core_circle(:,1)
>
> ;---Set points that are outside of the circle of data to missing
>
> lat_2d_outer_core:=lat_wrf(extrem_outer_core_inds
> (0,0):extrem_outer_core_inds(1,0),extrem_outer_core_inds(0,
> 1):extrem_outer_core_inds(1,1))
> lon_2d_outer_core:=lon_wrf(extrem_outer_core_inds
> (0,0):extrem_outer_core_inds(1,0),extrem_outer_core_inds(0,
> 1):extrem_outer_core_inds(1,1))
>
> in_circle_outer_core := gc_inout(lat_2d_outer_core,
> lon_2d_outer_core,clat_outer_core,clon_outer_core)
>
> size_in_circle_outer_core=dimsizes(in_circle_outer_core)
>
> ;;;;;;;;;;;;;;  on arc  ;;;;;;;;;;;
>
> clat_outer_core_size = dimsizes(clat_outer_core)
> print("clat_outer_core_size="+clat_outer_core_size)
>
> clat_outer_core_2d=reshape(clat_outer_core,(/floattointeger(clat_outer_
> core_size/2),2/))
> clon_outer_core_2d=reshape(clon_outer_core,(/floattointeger(clat_outer_
> core_size/2),2/))
>
> size_lat_2d_outer_core=dimsizes(lat_2d_outer_core)
>
> lat_2d_outer_core_reshape = reshape(lat_2d_outer_core,(/
> size_lat_2d_outer_core(0)*size_lat_2d_outer_core(1),2/))
> lon_2d_outer_core_reshape = reshape(lon_2d_outer_core,(/
> size_lat_2d_outer_core(0)*size_lat_2d_outer_core(1),2/))
>
> lat_2d_outer_core_to_1d = ndtooned(lat_2d_outer_core)
> lon_2d_outer_core_to_1d = ndtooned(lon_2d_outer_core)
>
> onarc_outer_core = gc_onarc(lat_2d_outer_core_to_
> 1d,lon_2d_outer_core_to_1d,clat_outer_core_2d,clon_outer_core_2d)
>
> ;;;;;;;;;;;;; on arc ;;;;;;;;;;;
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;;;;;;;;;;;;;;;;;;;
> ;;;;;;;;;   /\ /\ /\ /\ end of define circle sized structure /\ /\ /\ /\
> /\;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;;;;;;;;;;;;;;;;;;;
> end
>
>
>
>
> Best Regards
> =================================
> Nafiseh Pegahfar
> Assistant Professor
> Iranian National Institute for Oceanography and Atmospheric Science
> (http://www.inio.ac.ir)
> Phone: (0098)21- 66944873-5 Ext. 224
> Fax: (0098)21- 66944869
> Email: (pegahfar at inio.ac.ir)
> (pegahfar at ut.ac.ir)
> =================================
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180119/2b117cf0/attachment.html>


More information about the ncl-talk mailing list