Search code examples
rplotrasternetcdf

NetCDF display of trajectory data in R


I am trying to display the following classic NetCDF data Data Link Click Here. This data has 38 variables trajectory data and i am trying to extract one of them "air temp". Though i am able to extract the data but unfortunately i am not able to display the data in R using plot. The display should be like the following figure enter image description here, which i plotted in panoply. Another problem is that i am not able to stack all the data since the trajectories are all different. Is there any way i could display and stack all the files of one variable.

> library(ncdf4)
> library(rasterVis)
> library(raster)
lon <- ncvar_get(ncin, "lon")
lat <- ncvar_get(ncin, "lat")
data <-lon <- ncvar_get(ncin, "air_temp_ac")          #to extract variable
> dim(data)
[1] 6639
> dput(data)
structure(c(NA, 15, 14, 13, 11, NA, 11, 11, 11, 11, 11, NA, 14, 
14, 12, NA, 12, 14, 14, 14, 16, 21, 25, 27, 26, 19, 21, 22, 22, 
21, 22, 22, 23, 24, 22, 23, 25, 23, 26, 25, 25, 23, 24, 24, 25, 
28, 28, 29, 29, 28, 28, 26, 27, 27, 29, 31, 31, 34, 36, 37, 38, 
41, 37, 37, 39, 36, 29, 33, 38, 35, 36, 36, 36, 37, 37, 36, 37, 
34, 33, 34, 38, 37, 37, 37, 37, 37, 39, 39, 39, 40, 39, 39, 40,
> plot(lon,lat,data)
Error in plot.xy(xy, type, ...) : invalid plot type

I tried

> r <- raster(t(data), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
> r <- flip(r, direction='y')
> plot(r)

But, the plot is out of bound enter image description here


Solution

  • Here's a solution using ggplot2, sf and rnaturalearth. Using your code I couldn't extract the variable "air_temp_ac" from the nc file. The actual variable is "air_temp_AC" (note this is case sensitive)

    library(ncdf4)
    library(ggplot2)
    library(rnaturalearth) 
    library(sf)
    #> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
    
    world <- ne_countries(scale = "medium", returnclass = "sf")
    
    ncin  <- nc_open("~/nc.nc")
    
    df    <- data.frame(lon    = ncvar_get(ncin, "lon"), 
                        lat    = ncvar_get(ncin, "lat"),
                        Kelvin = ncvar_get(ncin, "air_temp_AC"))
    
    ggplot(data = world) +
      geom_rect(data = NULL, aes(xmin = 0, xmax= 90, ymin = 0, ymax = 60), fill = "#454565") +
      geom_sf(fill = "#8f9f7d") +
      coord_sf(xlim = c(0, 90), ylim = c(0, 55), expand = FALSE) +
      geom_point(data = df, aes(x = lon, y = lat, colour = Kelvin)) +
      scale_color_gradientn( colours = c("skyblue", "blue", "red", "orange", "yellow")) +
      labs(title = "Air temperature (K)", x = "longitude", y = "latitude")
    

    enter image description here

    Created on 2020-02-23 by the reprex package (v0.3.0)