How do I map NetCDF files in R?

General
, , ,
1

I am trying to plot a netCDF (.nc) file in R using ggplot to create a map. The data I'm using can be found (ESGF MetaGrid with variable "pr", frequency "mon" and source ID "EC-Earth3" File --> pr_Amon_EC-Earth3_historical_r1i1p1f1_gr_199801-199812)

I want to find out how to zoom in and add a title? Please help! I am so close!

Can I manipulate this to put it in ggplot to give me more tools?

Essentially, the goal is: to take a slice of the data (such as time = 1 is January) and have a U.S. and a global map of the differences in precipitation for that dataset.

I followed this tutorial: (RPubs - How to open and work with NetCDF data in R)

#https://rpubs.com/boyerag/297592

library(ncdf4) # package for netcdf manipulation
library(raster) # package for raster manipulation
library(rgdal) # package for geospatial analysis
library(ggplot2) # package for plotting

nc_data <- nc_open('pr_Amon_EC-Earth3_historical_r1i1p1f1_gr_199801-199812.nc')

# Save the print(nc) dump to a text file
{
  sink('pr_Amon_EC-Earth3_historical_r1i1p1f1_gr_199801-199812.txt')
  print(nc_data)
  sink()
}

lon <- ncvar_get(nc_data, "lon")
lon[lon > 180] <- lon[lon > 180] - 360


lat <- ncvar_get(nc_data, "lat", verbose = F)
t <- ncvar_get(nc_data, "time")

head(lon) 

ndvi.array <- ncvar_get(nc_data, "pr") # store the data in a 3-dimensional array
dim(ndvi.array) 

fillvalue <- ncatt_get(nc_data, "pr", "_FillValue")
fillvalue

nc_close(nc_data) 

ndvi.array[ndvi.array == fillvalue$value] <- NA
ndvi.slice <- ndvi.array[, , 1] 
dim(ndvi.slice)

r <- raster(t(ndvi.slice), 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)

#add countries!!!!

library(maptools)
data(wrld_simpl)
plot(wrld_simpl, add = TRUE)

#showing precip in january 1998

precip1998
precip1998605×538 37.9 KB

2

I don't know if it will help but you should be able to read netcdf file using sf package. This will ease data manipulation. ggplot2 offers a geom_sf function that could help too.
Otherwise, every spatial graph package would work I guess.

1 Like
3

thank you! do you have an example of how to read and manipulate netcdf files with the sf packages? I tried to look it up but wasn't able to find much.

5

Check out the stars package for manipulating spatiotemporal arrays. Within stars is the geom_stars function that allows for rasters to play nice with ggplot2.

Once you have a stars object in your environment, you can add layers to your ggplot in the same way you'd work with non-spatial data. Importantly, you'll need to make sure the layers have the same coordinate reference system.

In pseudocode (I can't find your data set):

library(sf)
library(ggplot2)
library(stars)
library(dplyr)

stars_object <- raster::raster("your_data_set.nc") %>% st_as_stars()
sf_object <- sf::st_read("second_data_set.shp")

#Make sure they have the same CRS
sf::st_crs(stars_object) <- sf::st_crs(sf_object)

#Crop to bounding box (insert coordinates in ...)
bb <- sf::st_bbox(xmin = ..., xmax = ..., ymin = ..., ymax = ...)
stars_object <- stars_object[bb]

#plot
ggplot()+
    geom_stars(data = stars_object) +
    geom_sf(data = sf_object)

Subsetting stars objects by layer is pretty straightforward if you follow the vignettes starting here.

Meant to reply to @KarmicDreamwork

3 Likes
6

Thanks! Any netcdf file will work from cmip6, they are all 4 dimensional.
I'm trying your code but I don't know how to put in the shapefile for the sf_object.

I downloaded a random world countries shapefile from Esri (https://www.arcgis.com/home/item.html?id=2ca75003ef9d477fb22db19832c9554f) and am met with this error:

library(sf)
library(ggplot2)
library(stars)
library(dplyr)

stars_object <- raster::raster(file.choose(), verbose = TRUE, write = FALSE) %>% st_as_stars()
sf_object <- sf::st_read("countries.shp")

> Cannot open data source countries.shp
> Error in CPL_read_ogr(dsn, layer, query, as.character(options), quiet,  : 
>   Open failed.

What am I doing wrong?

7

I'm not quite sure - I'm able to download and load the countries shapefile into my environment. Are you sure that you're specifying the correct (and full) file path?

8

Hmm, I was able to get the shapefile but when I map it, it looks like this. Do you know what's going on? Also, is there a way to just have the lines for the countries?

image
image1200×900 89 KB

I am getting an error when trying to set the boundaries as well:

#Crop to bounding box (insert coordinates in ...)
bb <- sf::st_bbox(xmin = -171.791110603, xmax = -66.96466, ymin = 18.91619, ymax = 71.3577635769)
stars_object <- stars_object[bb]
Error in is.numeric(bb) && length(bb) == 4 : 
  argument "obj" is missing, with no default

My code:

install.packages("stars")
library(sf)
library(ggplot2)
library(stars)
library(dplyr)

stars_object <- raster::raster("pr_Amon_EC-Earth3_historical_r1i1p1f1_gr_199801-199812.nc") %>% st_as_stars()
sf_object <- sf::st_read("countries.shp")

#Make sure they have the same CRS
sf::st_crs(stars_object) <- sf::st_crs(sf_object)

#U.S. 'US': ('United States', (-171.791110603, 18.91619, -66.96466, 71.3577635769))
#Crop to bounding box (insert coordinates in ...)
bb <- sf::st_bbox(xmin = -171.791110603, xmax = -66.96466, ymin = 18.91619, ymax = 71.3577635769)
stars_object <- stars_object[bb]

#plot
ggplot()+
  geom_stars(data = stars_object) +
  geom_sf(data = sf_object)
9

Also, for my map here, it looks like my data is set at 0-360 while the map is -180-180. anyone know how to fix this?

10

You can use the raster::rotate() function to convert your data from the 0-360 to -180-180 longitude scale prior to turning it into a stars object. You should then be able to plot the two data sets on the same scale via ggplot2.

My guess is that the mismatch in longitude scaling is what's preventing you from cropping the data set in your previous comment. Your precipitation raster starts at 0 ° longitude on the -180-180 scale

11

thank you. what would that look like, though?

12

Check out the documentation by entering ?raster::rotate() to see examples of how to use the function. raster::rotate() just takes a raster as an argument, so in your case it'd probably look like this:

#load the raster
r <- raster::raster("your_file.nc")

#rotate the raster
r_rotated <- raster::rotate(r)

hope that helps!

13

I'll try that. it's giving me an error "Error: must request at least one colour from a hue palette" and I can't find any solutions on google. I tried giving ggplot different colors but none of them are working. Any idea?

14

It looks like this post addresses the error you're reporting. What variables are being mapped to color or fill in the ggplot function? They may be empty or full of NAs.

I can't help you without being able to reproduce your error - try making a reprex if you can. I suggest including a link to a dataset on 0-360 longitude scale that leads to the same error.

15

Thank you i will check it out. Just got a new computer so i apologize for the delay.

Do you know how to zoom in on an area for the netcdf map?

16

Okay, I am having issues getting this to work. I looked at the other forum but it's just saying that my data is now messed up and full of NAs

I'm using EC-Earth GCM (under experiment id) for precipitation ("pr" under variable) variable for monthly (under frequency), year 2050 from the cmip6 website but any netcdf file would work (https://esgf-node.llnl.gov/search/cmip6/)

I got the shape file from the ArcGIS website.

> library(ncdf4)
> library(sf)
> library(ggplot2)
> library(stars)
> library(dplyr)
> library(maps)
> library(sp)
> library(maptools)
> library(raster)
> 


> 
> #attempt to rotate from 0-360 to -180 to 180 lon

> r <- raster::raster(file.choose())
#open pr_Amon_EC-Earth3_ssp370_r1i1p1f1_gr_205001-205012.nc
> r1 <- raster::rotate(r)
> 


> print(r1)
class      : RasterLayer 
dimensions : 256, 512, 131072  (nrow, ncol, ncell)
resolution : 0.703125, 0.7016692  (x, y)
extent     : -180, 180, -89.81366, 89.81366  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 
source     : memory
names      : Precipitation 
values     : -1.745207e-25, 0.000460271  (min, max)
z-value    : 2050-01-16 

> 
> stars_object <- raster::raster(r1) %>% st_as_stars()

> sf_object <- sf::st_read(file.choose())
#open Countries_WGS84.shp

Reading layer `Countries_WGS84' from data source `/Users/8kk/Desktop/RStudio Projects/NetCDF Mapping/Countries_WGS84/Countries_WGS84.shp' using driver `ESRI Shapefile'
Simple feature collection with 251 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.6236
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

> 
> #Make sure they have the same CRS
> sf::st_crs(stars_object) <- sf::st_crs(sf_object)
> 
> ggplot()+
+   geom_stars(data = stars_object) +
+   geom_sf(data = sf_object) # + aes(color = values)

**Error: Must request at least one colour from a hue palette.**


> str (stars_object)
List of 1
 $ layer: logi [1:512, 1:256] NA NA NA NA NA NA ...
 - attr(*, "dimensions")=List of 2
  ..$ x:List of 7
  .. ..$ from  : num 1
  .. ..$ to    : int 512
  .. ..$ offset: num -180
  .. ..$ delta : num 0.703
  .. ..$ refsys: chr "+proj=longlat +datum=WGS84 +no_defs"
  .. ..$ point : logi NA
  .. ..$ values: NULL
  .. ..- attr(*, "class")= chr "dimension"
  ..$ y:List of 7
  .. ..$ from  : num 1
  .. ..$ to    : int 256
  .. ..$ offset: num 89.8
  .. ..$ delta : num -0.702
  .. ..$ refsys: chr "+proj=longlat +datum=WGS84 +no_defs"
  .. ..$ point : logi NA
  .. ..$ values: NULL
  .. ..- attr(*, "class")= chr "dimension"
  ..- attr(*, "raster")=List of 3
  .. ..$ affine     : num [1:2] 0 0
  .. ..$ dimensions : chr [1:2] "x" "y"
  .. ..$ curvilinear: logi FALSE
  .. ..- attr(*, "class")= chr "stars_raster"
  ..- attr(*, "class")= chr "dimensions"
 - attr(*, "class")= chr "stars"
17

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.