Get a variable at specific location

The NCL graphics package.

Get a variable at specific location

Postby leroygr » Fri Jul 22, 2011 4:49 am

Hi everybody,

I'm a new user of WRF and NCL and I have a question...

I want to have the values of wind, temperature, atmospheric pressure and air density at a specific location, in form of a time serie in netcdf format. I have a WRF output but I don't know how to use NCL to get these variables at the location I want (X:574565, Y: 5585549 in UTM31 WGS84 coordinate system) and for the time series of my simulation... I saw a lot of NCL scripts (wrf_user_ll_tij for example) but I'm a little lost :( .
Is it also possible to interpolate these at a specific height?

Somebody can help me with this issue?

Thank you very much.

G.L
leroygr
 
Posts: 28
Joined: Wed Jul 20, 2011 10:55 am

Re: Get a variable at specific location

Postby jimmyc » Sun Jul 24, 2011 7:23 am

All of this is straightforward.
In order to extract data at a point use the function you referred to, to convert from lat-lon to grid indice (closest grid point). Then you can extract the time series from the model at that specific point.

To write a netcdf file:http://www.ncl.ucar.edu/Applications/o-netcdf.shtml

Yes you can interpolate. The wrf tutorial with NCL is a good place to start:
http://www.mmm.ucar.edu/wrf/OnLineTutorial/Graphics/NCL/index.html
The views expressed in this message do not necessarily reflect those of NOAA or the National Weather Service or the University of Oklahoma.
James Correia, Jr
jimmyc
 
Posts: 519
Joined: Tue Apr 15, 2008 1:10 am

Re: Get a variable at specific location

Postby leroygr » Tue Jul 26, 2011 2:11 am

Thanks for your answer!

I succeed to get windspeed at a specific location. Now i want to write the values for each time in an ASCII file. I used the command asciiwrite and I get the values in a ASCII output file. I want to go further and put more informations in the output file (latitude, longitude and time). I tried to apply example 4 from asciiwrite ncl documentation but it doesn't work.

Somebody have an idea to help me?

Here is my NCL script:

Code: Select all
;************************************************************************************
;* NCL script to generate netCDF file containing wind speed at specific location    *
;* Example for wrfout_d01_2000-01-24_12:00:00.nc                                    *
;*                                        *
;************************************************************************************


; --------------  LOAD FUNCTIONS AND PROCEDURES ----------------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

; --------------  BEGINING OF NCL SCRIPT ----------------


begin
;********************************************************
; read in netCDF file and make a loop for all time steps
;********************************************************
  in     = addfile("./wrfout_d01_2000-01-24_12:00:00.nc","r")


times = wrf_user_list_times(in)            ; get times in the file
ntimes = dimsizes(times)            ; number of times in the file
wind_speed = new(ntimes,float)            ; creation of a windspeed vector at each time step
;print(ntimes)

do it = 0,ntimes-1                  ;Loop for the time: it= starting time
print("Working on time " + it )
time = it


;************************************************
;  - Select lon & lat of the point of interest -
;************************************************


 res = True      
 res@returnInt = True               ; False : return real values, True: return interger values
 lat = 32.5                  ; Latitude of the point of interest
 lon = -87                  ; Longitude of the point of interest
 point = wrf_user_ll_to_ij(in,lon,lat,res)       ; wrf_user_ll_to_ij(nc_file,lon,lat,opt)
 
 x = point(1)
 y = point(0)

 print("X location is: " + x)            ; print the value of X at the screen
 print("Y location is: " + y)            ; print the value of Y at the screen


;************************************************
;  - extract wind components and unstagger them -
;************************************************


; ----- Hence U and V are staggered, we need to  average them to get the point on the mass grid using wrf_user_unstagger

  U     = wrf_user_getvar(in, "U", time)
  ua_in = wrf_user_unstagger(U,U@stagger)
  ua    = ua_in(0:0,y,x)             ;ua(bottom_up,grid_lon,grid_lat)
  V     = wrf_user_getvar(in, "V", time)
  va_in = wrf_user_unstagger(V,V@stagger)      
  va    = va_in(0:0,x,y)            ;va(bottom_up,grid_lon,grid_lat)
  wind_speed(it) = sqrt(ua^2+va^2)
  copy_VarCoords(ua,wind_speed(it))                         ; copy coord vars to speed





end do                     ; end of time loop
;************************************************
;  - Write wind speed in ascii file -
;************************************************




  XLAT = lat              ; (nlat)
  YLON = lon              ; (mlon)
  windspeed   = wind_speed             ; (ntim,nlat,mlon)

  dimwindspeed  = dimsizes(windspeed)
  ntim  = dimwindspeed(0)
  nXLAT  = dimwindspeed(1)
  mYLON  = dimwindspeed(2)

  npts  = nXLAT*mYLON        ; total number of grid points

  fName = "foo.txt"
  data  = new( npts, "string")

  npt   = -1

  do nl=0,nXLAT-1         
    do ml=0,mYLON-1

       npt  = npt + 1   
       data(npt) = sprinti("%0.5i", (npt+1) ) 
       data(npt) = data(npt) + sprintf("%7.1f ",XLAT(nl))
       data(npt) = data(npt) + sprintf("%7.1f ",YLON(ml))

      do nt=0,ntim-1
         data(npt) = data(npt) + sprintf("%10.3f ", windspeed(nt,nl,ml))
      end do
    end do
  end do

  asciiwrite (fName , data)


; --------------  END OF NCL SCRIPT ----------------

end


Thank you very much!

Greg
leroygr
 
Posts: 28
Joined: Wed Jul 20, 2011 10:55 am

Re: Get a variable at specific location

Postby jimmyc » Wed Jul 27, 2011 5:35 pm

what was the error?
The views expressed in this message do not necessarily reflect those of NOAA or the National Weather Service or the University of Oklahoma.
James Correia, Jr
jimmyc
 
Posts: 519
Joined: Tue Apr 15, 2008 1:10 am

Re: Get a variable at specific location

Postby leroygr » Thu Jul 28, 2011 1:47 am

Here is what i get when i run the script:

Code: Select all
[greg@Uranus:~/WRF/post_processing/NCL/scripts]$ ncl wind_speed2.ncl
 Copyright (C) 1995-2011 - All Rights Reserved
 University Corporation for Atmospheric Research
 NCAR Command Language Version 6.0.0
 The use of this software is governed by a License Agreement.
 See http://www.ncl.ucar.edu/ for more details.


Variable: times
Type: string
Total Size: 20 bytes
            5 values
Number of Dimensions: 1
Dimensions and sizes:   [5]
Coordinates:
Number Of Attributes: 2
  description :   times in file
  _FillValue :   missing
(0)   2000-01-24_12:00:00
(1)   2000-01-24_15:00:00
(2)   2000-01-24_18:00:00
(3)   2000-01-24_21:00:00
(4)   2000-01-25_00:00:00
(0)   Working on time 0
(0)   X location is: 19
(0)   Y location is: 21
(0)   Working on time 1
(0)   X location is: 19
(0)   Y location is: 21
(0)   Working on time 2
(0)   X location is: 19
(0)   Y location is: 21
(0)   Working on time 3
(0)   X location is: 19
(0)   Y location is: 21
(0)   Working on time 4
(0)   X location is: 19
(0)   Y location is: 21
fatal:Subscript out of range, error in subscript #0
fatal:An error occurred reading dimwindspeed
fatal:Execute: Error occurred at or near line 87 in file wind_speed2.ncl


Thank you!
leroygr
 
Posts: 28
Joined: Wed Jul 20, 2011 10:55 am

Re: Get a variable at specific location

Postby mallesh4science » Tue Jan 06, 2015 8:23 am

Dear Leroygr,

I am also facing the same problem. How did you resolve it.

Regards,
Malleswararao
mallesh4science
 
Posts: 11
Joined: Wed May 22, 2013 4:49 am

Re: Get a variable at specific location

Postby bradyhalvorson » Thu Jun 21, 2018 9:21 am

I'm not sure if you ever solved this for yourself but I figured I would post my solution to help others.

I believe the error was due to the wind_speed variable only having 1 dimension (the data itself) and not 3 as was assumed in your code. If you set nXLAT = 1 and mYLON = 1 and then change line 106 to have nt as the only input in windspeed, the code runs fine for me. This seemed like a simplification but you're only interested in data at one point so I'm not sure why you would need to carry over the location data from U because the location is what you specified at the beginning.

I also had to change your variable wind_speed to wind_speed_var due to the fact that wind_speed is now a function in NCL and will produce an error otherwise.

Code: Select all
;************************************************************************************
;* NCL script to generate netCDF file containing wind speed at specific location    *
;* Example for wrfout_d01_2000-01-24_12:00:00.nc                                    *
;*                                        *
;************************************************************************************


; --------------  LOAD FUNCTIONS AND PROCEDURES ----------------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

; --------------  BEGINING OF NCL SCRIPT ----------------


begin
;********************************************************
; read in netCDF file and make a loop for all time steps
;********************************************************
  in     = addfile("./wrfout_d01_2018-06-19_06:00:00.nc","r")


times = wrf_user_list_times(in)            ; get times in the file
ntimes = dimsizes(times)            ; number of times in the file
wind_speed_var = new(ntimes,float)            ; creation of a windspeed vector at each time step
;print(ntimes)

do it = 0,ntimes-1                  ;Loop for the time: it= starting time
print("Working on time " + it )
time = it


;************************************************
;  - Select lon & lat of the point of interest -
;************************************************


 res = True
 res@returnInt = True               ; False : return real values, True: return interger values
 lat = 44.3                  ; Latitude of the point of interest
 lon = -94.4                 ; Longitude of the point of interest
 point = wrf_user_ll_to_ij(in,lon,lat,res)       ; wrf_user_ll_to_ij(nc_file,lon,lat,opt)

 x = point(1)
 y = point(0)


 print("X location is: " + x)            ; print the value of X at the screen
 print("Y location is: " + y)            ; print the value of Y at the screen


;************************************************
;  - extract wind components and unstagger them -
;************************************************


; ----- Hence U and V are staggered, we need to  average them to get the point on the mass grid using wrf_user_unstagger

  U     = wrf_user_getvar(in, "U", time)
  ua_in = wrf_user_unstagger(U,U@stagger)
  ua    = ua_in(0:0,y,x)             ;ua(bottom_up,grid_lon,grid_lat)
  V     = wrf_user_getvar(in, "V", time)
  va_in = wrf_user_unstagger(V,V@stagger)
  va    = va_in(0:0,x,y)            ;va(bottom_up,grid_lon,grid_lat)
  wind_speed_var(it) = sqrt(ua^2+va^2)
  copy_VarCoords(ua,wind_speed_var(it))                         ; copy coord vars to speed





end do                     ; end of time loop
;************************************************
;  - Write wind speed in ascii file -
;************************************************




  XLAT = lat              ; (nlat)
  YLON = lon              ; (mlon)
  windspeed   = wind_speed_var             ; (ntim,nlat,mlon)

  dimwindspeed  = dimsizes(windspeed)
  ntim  = dimwindspeed(0)
  nXLAT  = 1
  mYLON  = 1

  npts  = nXLAT*mYLON        ; total number of grid points

  fName = "foo.txt"
  data  = new( npts, "string")

  npt   = -1

  do nl=0,nXLAT-1
    do ml=0,mYLON-1

       npt  = npt + 1
       data(npt) = sprinti("%0.5i", (npt+1) )
       data(npt) = data(npt) + sprintf("%7.1f ",XLAT(nl))
       data(npt) = data(npt) + sprintf("%7.1f ",YLON(ml))

      do nt=0,ntim-1
         data(npt) = data(npt) + sprintf("%11.3f ", windspeed(nt))
      end do
    end do
  end do

  asciiwrite (fName , data)


; --------------  END OF NCL SCRIPT ----------------

end
bradyhalvorson
 
Posts: 1
Joined: Thu Jun 21, 2018 8:57 am


Return to NCL

Who is online

Users browsing this forum: No registered users and 1 guest