[ncl-talk] Extracting A Time Series of Forcing Data Information from Specific Locations

Smith, Stuart smit1770 at purdue.edu
Fri Mar 9 15:29:41 MST 2018


Not a problem. Thank you for your assistance Dennis, and providing feedback on my mistakes. This is very helpful. Thank you for your time, and have a nice day.

 

Regards,

 

-Stuart

 

From: Dennis Shea [mailto:shea at ucar.edu] 
Sent: Friday, March 9, 2018 4:40 PM
To: Smith, Stuart <smit1770 at purdue.edu>
Cc: Ncl-talk <ncl-talk at ucar.edu>
Subject: Re: [ncl-talk] Extracting A Time Series of Forcing Data Information from Specific Locations

 

Apologies: I was in too much of a rush.

The following uses the files you sent.
======================

; vector of names of all of the files and add full file path
  
diri = "./"
fili = systemfunc("ls "+diri+"NLDAS_FORA0125_H.A2009*.002.grb")
print(fili)   ; make sure the desired files are listed in the correct order

; read in all files

;ListSetType (grb_file, "join")   
setfileoption("grb","SingleElementDimensions","Initial_time") ; initial_time0_hours

grb_file = addfiles(fili,"r")
data = grb_file[:]->DSWRF_110_SFC
data!0 = "time"
printVarSummary(data)

;====
   LAT  = (/ 48.5625/)   ; read from text file
   LON  = (/-98.9375/)
   LAT at units = "degrees_north"
   LON at units = "degrees_east"

; Simple approx:

   print(data(:,{LAT},{LON}) )
;=====
; or better use:

;https://www.ncl.ucar.edu/Document/Functions/Contributed/linint2_points_Wrap.shtml

points = linint2_points(data&lon_110,data&lat_110,data, False, LON,LAT, 0)
printVarSummary(points)
print("--------------")

+++++++++++++++++++++++++++++
Variable: fili
Type: string
Total Size: 24 bytes
            3 values
Number of Dimensions: 1
Dimensions and sizes:    [3]
Coordinates: 
(0)    ./NLDAS_FORA0125_H.A20090316.0100.002.grb
(1)    ./NLDAS_FORA0125_H.A20090316.0200.002.grb
(2)    ./NLDAS_FORA0125_H.A20090316.0300.002.grb

Variable: data
Type: float
Total Size: 1247232 bytes
            311808 values
Number of Dimensions: 3
Dimensions and sizes:    [time | 3] x [lat_110 | 224] x [lon_110 | 464]
Coordinates: 
            time: [1833841..1833843]
            lat_110: [25.063..52.938]
            lon_110: [-124.938..-67.063]
Number Of Attributes: 12
  sub_center :    NESDIS Office of Research and Applications
  center :    US National Weather Service - NCEP (WMC)
  long_name :    Downward shortwave radiation flux
  units :    W/m^2
  _FillValue :    1e+20
  level_indicator :    1
  gds_grid_type :    0
  parameter_table_version :    130
  parameter_number :    204
  model :    MESO NAM Model (currently 12 km)
  forecast_time :    0
  forecast_time_units :    hours
------
Variable: data (subsection)       <====  print(data(:,{LAT},{LON}) )
Type: float
Total Size: 12 bytes
            3 values
Number of Dimensions: 1
Dimensions and sizes:    [time | 3]
Coordinates: 
            time: [1833841..1833843]
Number Of Attributes: 14
  lon_110 :    -98.938
  lat_110 :    48.563
  forecast_time_units :    hours
  forecast_time :    0
  model :    MESO NAM Model (currently 12 km)
  parameter_number :    204
  parameter_table_version :    130
  gds_grid_type :    0
  level_indicator :    1
  _FillValue :    1e+20
  units :    W/m^2
  long_name :    Downward shortwave radiation flux
  center :    US National Weather Service - NCEP (WMC)
  sub_center :    NESDIS Office of Research and Applications
==========

Variable: points
Type: float
Total Size: 12 bytes
            3 values
Number of Dimensions: 2
Dimensions and sizes:    [time | 3] x [pts | 1]
Coordinates: 
            time: [1833841..1833843]
            pts: [0..0]
Number Of Attributes: 14
  _FillValue :    1e+20
  sub_center :    NESDIS Office of Research and Applications
  center :    US National Weather Service - NCEP (WMC)
  long_name :    Downward shortwave radiation flux
  units :    W/m^2
  level_indicator :    1
  gds_grid_type :    0
  parameter_table_version :    130
  parameter_number :    204
  model :    MESO NAM Model (currently 12 km)
  forecast_time :    0
  forecast_time_units :    hours
  xcoord :    -98.9375
  ycoord :    48.5625








 

 

On Fri, Mar 9, 2018 at 1:59 PM, Dennis Shea <shea at ucar.edu <mailto:shea at ucar.edu> > wrote:

[1]

The following is incorrect. grb_list is a one-dimensional *array* of file names containing strings (text).  You have it as a single string.

 

    grb_file = addfiles("grb_list","r")  

 

Use 

 

    grb_file = addfiles( grb_list ,"r")  ; pass the 1-d array of file names

=============
[2]

You want time series, N'est ce pas? Remove the following:

    ListSetType (grb_file, "join")            ; NO: NOT if you want a time serie
=============
[3]


diri = ncargpath("/depot/phig/apps/LIS_VIC/NLDAS_Forcing/discover/nobackup/projects/lis/MET_FORCING/NLDAS2.FORCING/2009/075/")

grb_file = addfiles("grb_list","r")

data      = grb_file[:]->DSWRF_110_SFC

printVarSummary(data)

 

ariable: data
Type: float
Total Size: 2494464 bytes
            623616 values
Number of Dimensions: 2
Dimensions and sizes:    [lat_110 | 1344] x [lon_110 | 464]
Coordinates: 
            lat_110: [25.063..52.938]
            lon_110: [-124.938..-67.063]

 

=====

   LAT  = (/ 48.5625 /)   ; read from text file

   LAT at units = "degrees_north"

   LON = (/-98.9375/)

   LON&units = "degrees_east"

   npts = dimsizes(LAT)

 

; Simple approx:

 

   print( data(:,{LAT),{LON}) )

=====

 

or better use:

https://www.ncl.ucar.edu/Document/Functions/Contributed/linint2_points_Wrap.shtml

 

points = linint2_points(data&lon,data&lat,data, False, LON,LAT, 0)
printVarSummary(points)
print("--------------")

 

Good Luck

 

On Fri, Mar 9, 2018 at 9:37 AM, Smith, Stuart <smit1770 at purdue.edu <mailto:smit1770 at purdue.edu> > wrote:

Hello,

 

I would like to compare two sets of forcing data. One is in a .grb file while the other is a .txt. When working the GRIB files (1 hour timesteps), I would like my .ncl script to reads in a time series of data, with the goal to plot specific weather data variables from specific locations to compare the .txt files. 

 

However, when running my .ncl script listed below I receive a Segmentation fault error. I’m not sure if this is because I am not defining the latitude and longitude coordinates soon enough, causing an error by reading in too much data. I do not know how to print or plot specific latitude (48.5625) and longitude (-98.9375) points and print their forcing variables. Could you please help improve my .ncl script? I have also uploaded an example of the .grb files I am working with. Thank you for your time.

 

 

begin

 

;Go to directory of interest 

 

diri = ncargpath("/depot/phig/apps/LIS_VIC/NLDAS_Forcing/discover/nobackup/projects/lis/MET_FORCING/NLDAS2.FORCING/2009/075/")

 

 

;list names of all of the files and add ful file path

grb_list = systemfunc("ls "+diri+"NLDAS_FORA0125_H.*.002.grb") 

  

;print(grb_list)   ; make sure the desired files are listed in the correct order

 

; read in all files

 

grb_file = addfiles("grb_list","r")

ListSetType (grb_file, "join")

data = grb_file[:]->DSWRF_110_SFC

data at lat=grb_file[:]->lat_110  ;Lat. of interest 48.5625

data at lon=grb_file[:]->lon_110  ;Lon. of interest -98.9375

 

print(data)

 

 

end

 

 

Regards,

 

-Stuart 

 

 

_______________________________________________
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

 

 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180309/d432f60f/attachment.html>


More information about the ncl-talk mailing list