[ncl-talk] Extracting Output Results from Specific Grid Cells

[Mary Haley]

Stuart,

Sorry for the delay.

The answer follows, but I want to point out that almost half of your
lat/lon values on the file are missing. I'll discuss this in a little bit.

If you look at one of your data files using ncl_filedump, you will see that
the lat/lon arrays are two-dimensional (I edited the output for clarity):

    Variable: lat
    Type: float
    Number of Dimensions: 2
    Dimensions and sizes: [ 29 <north_south> x 45 <east_west> ]
    Coordinates:
        Number of Attributes:        9
            units : degree_north
            standard_name : latitude
            long_name : latitude


    Variable: lon
    Type: float
    Number of Dimensions: 2
    Dimensions and sizes: [ 29 <north_south> x 45 <east_west> ]
    Coordinates:
        Number of Attributes:        9
            units : degree_east
            standard_name : longitude
            long_name : longitude

Data that is defined by two-dimensional lat/lon arrays is said to be on a
"curvilinear" grid. In order to sub-select regions from curvilinear data,
you need to use either getind_latlon2d or region_ind.

Also, you don't want to use the [:] syntax to read the lat/lon arrays off
the file, because then you will get back 2 x 29 x 45 arrays, which is
unnecessary if your lat/lon arrays don't change with each time step.

Instead, simply read off one timestep of the lat / lon using [0]:


  a         = addfiles(fils,"r")
  ListSetType (a, "join")
  data      = a[:]->Tair_f_tavg
  data at time = a[:]->time
  lon2d     = a[0]->lon
  lat2d     = a[0]->lat

You can then get the index locations into lat2d and lon2d that are equal to
LAT and LON as follows:

  LAT  = (/ 48.5625/)
  LON  = (/-98.9375/)

  nm = getind_latlon2d (lat2d,lon2d, LAT, LON)
  n = nm(0,0)
  m = nm(0,1)
  print("==================================================")
  print("Original lat/lon : " + LAT       +"   "+LON)
  print("Actual lat/lon   : " + lat2d(n,m)+"   "+lon2d(n,m))
  print("==================================================")

To slice the data at that particular location:


  data_slice = data(:,n,m)
  printVarSummary(data)
  printVarSummary(data_slice)
  printMinMax(data_slice,0)


We have some examples of using these functions. Please see:

http://www.ncl.ucar.edu/Applications/latlon_subset.shtml

These examples show the difference between region_ind and getind_latlon2d,
which is important if you want to grab a sub-region of data, rather than
just a single slice.

Back to the lat/lon values being missing. It's always incredibly important
to look at your data and understand it before you start doing calculations
on it, or try to plot it.  Dennis Shea mentioned to me offline that your
lat/lon arrays had a high number of missing values, so I simply added the
following two lines to check this:

  print("==================================================")
  print("# of valid latitude values    = "  + num(.not.ismissing(lat2d)))
  print("# of missing latitude values  = "  + num(ismissing(lat2d)))
  print("==================================================")
  print("# of valid longitude values   = " + num(.not.ismissing(lon2d)))
  print("# of missing longitude values = " + num(ismissing(lon2d)))
  print("==================================================")

which gave me the following output:

==================================================
# of valid latitude values    = 703
# of missing latitude values  = 602
==================================================
# of valid longitude values   = 703
# of missing longitude values = 602
==================================================

It's likely that the lat/lon arrays have missing values in the same
location.

Now maybe this is perfectly valid for your data, but it's something we
wanted to be sure you were aware of. It's somewhat unusual to see this many
missing values in lat / lon arrays.

--Mary





On Sat, Apr 28, 2018 at 12:32 PM, Smith, Stuart <smit1770 at purdue.edu> wrote:

> Hello,
>
>
>
> Attached is a NCL script I have been using to output model simulation
> results.  My goal is to enhance this script to plot outputs from specific
> grid cells.
>
>
>
> The plotting script was originally designed to read in netCDF files, which
> were outputted from the model simulation in a daily time step as a netCDF
> file  (ex. 200903180000.d01.nc). In order to create a timeseries, the
> script joins the files into an array and plots the total average data. I
> would like to enhance the NCL script to plot a time series of the raw
> outputted data from a defined grid cell.
>
>
>
> I have been able to do this successfully looking at individual files, but
> when I try joining the data into an array for time series analysis of the
> grid cell I get the following error “fatal:Dimension (east_west) of (data)
> does not have an associated coordinate variable fatal" on line 28. Could I
> please have assistance on understanding the error and correcting this
> issue? Attached are example output results from a model simulation as well.
> Thank you for your time.
>
>
>
> Regards,
>
>
>
> -Stuart
>
> _______________________________________________
> 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/20180502/d0b5eff2/attachment.html>