Search code examples
rggplot2rasternetcdfr-raster

Set Geographic Coordinate System to WGS84 for netCDF Raster Data in R

I have a netCDF file containing bathymetry data that I am trying to combine on the same plot with shapefiles that are in WGS84. When I plot the netCDF data it doesn't have any coordinate system, how can I set a coordinate system for the netCDF data so I can plot all my spatial files together on the same plot?

library(ncdf4) 
library(raster) 
library(rgdal) 
library(ggplot2)

# READ IN AND PARSE BATHYMETRY DATA
# from https://www.ngdc.noaa.gov/thredds/catalog/regional/catalog.htmldataset=regionalDatasetScan
#/albemarle_sound_s010_30m.nc

nc_data <- nc_open('albemarle_sound_s010_30m.nc')

# Get Data from netCDF File
x_coord <- ncvar_get(nc_data, "x")
y_coord <- ncvar_get(nc_data, "y")
Band1 <- ncvar_get(nc_data, "Band1")

# Close connection
nc_close('albemarle_sound_s010_30m.nc')

# Melt data into dataframe for plotting in ggplot2
Band1_df = melt(Band1)

# Plot Data
ggplot() + 
  geom_raster(aes(x = X1, y = X2, fill = value), data = Band1_df) +
  scale_fill_continuous(na.value = "transparent") +
  ylab("") +
  xlab("") +
  theme(panel.grid.major = element_line(color = gray(.5), linetype = "dashed", size = 0.7),
    panel.background = element_rect(fill = "black"),
    plot.title = element_text(color = "black", size = 16, hjust = 0.5))
Solution

  • The URL you provide is not valid. But I think this is the file you are using

    furl <- "https://www.ngdc.noaa.gov/thredds/fileServer/regional/albemarle_sound_s010_30m.nc"
f = basename(furl)
download.file(furl, f, mode="wb")

    With a ncdf file you can do

    library(terra)
x <- rast(f)
x
#class       : SpatRaster 
#dimensions  : 3536, 4013, 1  (nrow, ncol, nlyr)
#resolution  : 30, 30  (x, y)
#extent      : 325249.1, 445639.1, 3945457, 4051537  (xmin, xmax, ymin, ymax)
#coord. ref. : NAD83 / UTM zone 18N (EPSG:26918) 
#source      : albemarle_sound_s010_30m.nc 
#varname     : Band1 (GDAL Band Number 1) 
#name        : Band1 
 
plot(x)

    enter image description here

    And now you could use lines() and points() to add spatial vector data ("shapefiles"). or use ggplot or a more mapping oriented variatoin of that instead such as tmap.

    As you can see above, the raster data have a "UTM zone 18N" coordinate reference system (CRS). You say that you want to use longitude/latitude (I assume that is what you mean with WGS84). You can do that, but it makes more sense to transform the vector data to the coordinate reference system of the raster data.