[ncl-talk] How to extract data just over the circle defined by "nggcog" function
Kevin Hallock
hallock at ucar.edu
Fri Mar 16 19:02:18 MDT 2018
Hi Nafiseh,
I’m sorry for the delay in responding to you.
I’ve looked into gc_onarc() to determine how it is supposed to be used, and in the interest of clarifying things for you and others I’d like to share my experience.
For the calculation you’re trying to do, the two rightmost arguments represent arc segments of the circle generated by nggcog(). For example, if clon_outer_core was (/0, 1, 2, 3/) after running nggcog(), your clon_outer_core_2d should be (/(/0, 1/), (/1, 2/), (/2, 3/), (/3, 0/)/); each of the four pairs of longitudes would represent the start- and end-points of the great circle arcs that approximate a circle around lat_location and lon_location. The method you used below (using reshape()) will actually skip half of the arcs and would look like this: (/(/0, 1/), (/2, 3/)/). The following snippet of code should generate the necessary clat_outer_core_2d and clon_outer_core_2d arrays containing all arcs:
clat_outer_core_2d = new((/clat_outer_core_size, 2/), typeof(clat_outer_core))
clon_outer_core_2d = new((/clon_outer_core_size, 2/), typeof(clon_outer_core))
clat_outer_core_2d(:,0) = clat_outer_core
clat_outer_core_2d(0:clat_outer_core_size-2,1) = clat_outer_core(1:)
clat_outer_core_2d(clat_outer_core_size-1,1) = clat_outer_core(0)
clon_outer_core_2d(:,0) = clon_outer_core
clon_outer_core_2d(0:clat_outer_core_size-2,1) = clon_outer_core(1:)
clon_outer_core_2d(clat_outer_core_size-1,1) = clon_outer_core(0)
The sizes of the arguments passed to gc_onarc() must match because gc_onarc checks to see if the nth point in the first two arguments is on the arc defined by the nth index of the last two arguments. This means that if you want to check if a particular point is on any of the 144 arcs, your first two arguments will need to be arrays containing the same value 144 times. Mary’s solution below almost works, but needs to be adjusted to account for having the same value 144 times:
npts = dimsizes(lat_2d_outer_core_to_1d)
onarc_outer_core = new(npts,logical)
do n=0,npts-1
conform_lat = conform_dims(clat_outer_core_size, lat_2d_outer_core_to_1d(n), -1)
conform_lon = conform_dims(clat_outer_core_size, lon_2d_outer_core_to_1d(n), -1)
onarc_outer_core(n) = any(gc_onarc(conform_lat,conform_lon,clat_outer_core_2d,clon_outer_core_2d))
end do
Note that I also added the “any” function around the gc_onarc call — this is because the output of gc_onarc is a logical/boolean array “dimsizes(clat_outer_core_2d)” long which tells you which arc index your point was on. I suspect in this case that you don’t care which particular arc segment the point was on, just if it was on any arc segment of the circle.
Also, I investigated the source code for gc_inout() and found that it treats any points that are “on arc” as being “in”; under the hood, gc_inout() calls the same Fortran function used by gc_onarc() in addition to its own “in/out” checks. As such, I believe you could actually eliminate the gc_onarc section of your program.
Finally, I have some ideas to make this program a little simpler and possibly quicker and more accurate. Because the shape you’re using for gc_inout() and gc_onarc() is a circle instead of an arbitrary polygon, instead of calling those functions you could simply see if the distance between two points is less than or equal to a specified radius in degrees:
distance = gc_latlon(lat_2d_outer_core_to_1d(n), lon_2d_outer_core_to_1d(n), lat_location, lon_location, 2, 2)
if(distance .le. radius_outer_core) then
if(distance .eq. radius_outer_core) then
print("on")
else
print("in")
end if
else
print("out")
end if
I have attached an NCL script that defines “on_circle” and “in_circle” functions that mimic the behavior of gc_onarc and gc_inout, respectively.
Please let me know if you have any further questions.
I hope this helps!
Kevin
> On Jan 19, 2018, at 2:35 PM, Mary Haley via ncl-talk <ncl-talk at ucar.edu> wrote:
>
> 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 <mailto: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 <http://www.inio.ac.ir/>)
> Phone: (0098)21- 66944873-5 Ext. 224
> Fax: (0098)21- 66944869
> Email: (pegahfar at inio.ac.ir <mailto:pegahfar at inio.ac.ir>)
> (pegahfar at ut.ac.ir <mailto:pegahfar at ut.ac.ir>)
> =================================
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk <http://mailman.ucar.edu/mailman/listinfo/ncl-talk>
>
>
> _______________________________________________
> 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/20180316/815dab22/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: on_circle_example.ncl
Type: application/octet-stream
Size: 3910 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180316/815dab22/attachment.obj>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180316/815dab22/attachment-0001.html>
More information about the ncl-talk
mailing list