0

I would like to use the rasterVis package to plot contours of spatial data (e.g. using levelplot as in this example). However, the dataset comes from a NetCDF file with an irregular grid, as below:

lon(y,x) lat(y,x) var(y,x)

and there is no implication to the projection implied by lon/lat.

Is there a method to plot directly the dataset as raster data as in these figures without interpolation?

The headers of raster data include an extension of the grid and projection specification, which does not fit with my problem. Two-dimensional lon/lat arrays are not recognized by raster as coordinate systems.

Code and plots as in this example but with netcdf file as:

float lat(y, x) ;
    lat:standard_name = "latitude" ;
    lat:long_name = "Latitude" ;
    lat:units = "degrees_north" ;
    lat:nav_model = "grid_T" ;
float lon(y, x) ;
    lon:standard_name = "longitude" ;
    lon:long_name = "Longitude" ;
    lon:units = "degrees_east" ;
    lon:nav_model = "grid_T" ;
float icethic(time_centered, y, x) ;
    icethic:standard_name = "sea_ice_thickness" ;
    icethic:long_name = "Ice thickness (cell average)" ;
    icethic:units = "m" ;
    icethic:online_operation = "average" ;
    icethic:_FillValue = 1.e+20f ;
    icethic:missing_value = 1.e+20f ;
    icethic:coordinates = "time_centered nav_lon nav_lat" ;
Majid
  • 1,836
  • 9
  • 19
Andrea
  • 1
  • 1
    Hi Andrea. Welcome to SO! You need to minimise your questions to one at the time for each post so you can get the best answer. Plus, try to clarify the question as much as possible by explaining what you did and what your expected output is. – Majid May 07 '19 at 07:59

1 Answers1

0

thanks for quick feedback. I would like to plot with rasterVis a NetCDF with irregular grid (lon, lat are 2d arrays): 

netcdf temp {                                                                                                                                             
dimensions:                                                                                                                                               
    y = 292 ;                                                                                                                                         
    x = 362 ;                                                                                                                                         
    time_counter = UNLIMITED ; // (1 currently)
variables:
    float lat(y, x) ;
            lat:standard_name = "latitude" ;
            lat:long_name = "Latitude" ;
            lat:units = "degrees_north" ;
            lat:_CoordinateAxisType = "Lat" ;
    float lon(y, x) ;
            lon:standard_name = "longitude" ;
            lon:long_name = "Longitude" ;
            lon:units = "degrees_east" ;
            lon:_CoordinateAxisType = "Lon" ;
    double time_counter(time_counter) ;
            time_counter:standard_name = "time" ;
            time_counter:units = "days since 0-00-00 00:00:00" ;
            time_counter:calendar = "proleptic_gregorian" ;
    float votemper(time_counter, y, x) ;
            votemper:standard_name = "Temperature" ;
            votemper:long_name = "Temperature" ;
            votemper:units = "C" ;
            votemper:coordinates = "lon lat time_counter" ;
            votemper:_FillValue = 9.96921e+36f ;
            votemper:missing_value = 9.96921e+36f ;
            votemper:online_operation = "ave(x)" ;
            votemper:interval_operation = 3600.f ;
            votemper:interval_write = 2678400.f ;
            votemper:offline_operation = "ave(x)" ;
}

The code inspired by the rasterVis guide looks like:

library(raster)
library(rasterVis)
stackSIS <- stack("temp.nc")
idx <- c(as.Date('2008-01-15'))
SISmm <- setZ(stackSIS, idx)
names(SISmm) <- month.abb[1]
SISmm
levelplot(SISmm)

but the plot does not consider the lon/lat geographical coordinates as axes but the x,y indices of the array. Indeed when I ask the summary of the raster object I get:

class       : RasterStack 
dimensions  : 292, 362, 105704, 1  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 0.5, 362.5, 0.5, 292.5  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
names       : Jan 
time        : 2008-01-15

i.e. the "extent" considers the indexes and not the coordinates.

Thanks

Andrea
  • 1