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

[nafiseh pegahfar]

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)
=================================
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180118/b4dd45f3/attachment.html>