0

Issue:

I have a code in R that extracts monthly sea surface temperature (SST) values from a single Aqua Modis netCDF file (see below). However, I have a batch of 59 Aqua Modis netCDF files in one folder.

Aim:

My aim is to extract the variable's longitude, latitude, and SST from every netCDF file across all 59 netCDF files, convert them into raster files using the function stack::raster(), and then process these files.

My data frame has 650 rows, which are dolphin IDs. I would like to extract the average SST for each dolphin ID over the time period of 2016-2021. Once I have extracted the average value SST per dolphin ID, I would then like to write these values into a .csv file named as a column called ' Average_SST'.

Unfortunately, I cannot share my data because of ownership issues.

I am a complete novice at spatial analysis and I have been trying to solve this conundrum for nearly 4 days by watching YouTube, and reading tutorials and Stack Overflow. I now feel really confused and I believe I will have to loop through all the files to achieve my objective.

If anyone can help, I would like to thank you in advance.

R-code:

Open all 59 Aqua Modis netCDF files from their folder

##Open packages needed for our analysis
library(ncdf4)
library(terra)
library('RNetCDF')
library(raster)
library(sp)

filenames = list.files('~/Documents/Ocean_ColorSST_2016_2021',pattern='*.nc',full.names=TRUE)

##Open the 70 Aqua Modis netCDF files from their folder
SST <- nc_open(filenames[59])

##Print the results for SST
print(SST)

Results:

     3 variables (excluding dimension variables):
        short sst[lon,lat]   (Chunking: [87,44])  (Compression: shuffle,level 4)
            long_name: Sea Surface Temperature
            scale_factor: 0.00499999988824129
            add_offset: 0
            units: degree_C
            standard_name: sea_surface_temperature
            _FillValue: -32767
            valid_min: -1000
            valid_max: 10000
            display_scale: linear
            display_min: -2
            display_max: 45
        unsigned byte qual_sst[lon,lat]   (Chunking: [87,44])  (Compression: shuffle,level 4)
            long_name: Quality Levels, Sea Surface Temperature
            _FillValue: 255
            valid_min: 0
            valid_max: 5
        unsigned byte palette[eightbitcolor,rgb]   (Contiguous storage)  

     4 dimensions:
        lat  Size:4320 
            long_name: Latitude
            units: degrees_north
            standard_name: latitude
            _FillValue: -999
            valid_min: -90
            valid_max: 90
        lon  Size:8640 
            long_name: Longitude
            units: degrees_east
            standard_name: longitude
            _FillValue: -999
            valid_min: -180
            valid_max: 180
        rgb  Size:3 (no dimvar)
        eightbitcolor  Size:256 (no dimvar)
[1] ">>>> WARNING <<<  attribute data_bins is an 8-byte value, but R"
[1] "does not support this data type. I am returning a double precision"
[1] "floating point, but you must be aware that this could lose precision!"

    61 global attributes:
        product_name: AQUA_MODIS.20210901_20210930.L3m.MO.SST.sst.4km.nc
        instrument: MODIS
        title: MODISA Level-3 Standard Mapped Image
        project: Ocean Biology Processing Group (NASA/GSFC/OBPG)
        platform: Aqua
        temporal_range: month
        processing_version: R2019.0
        date_created: 2021-12-03T08:21:22.000Z
        history: l3mapgen par=AQUA_MODIS.20210901_20210930.L3m.MO.SST.sst.4km.nc.param 
        l2_flag_names: LAND,HISOLZEN
        time_coverage_start: 2021-09-01T00:45:00.000Z
        time_coverage_end: 2021-10-01T02:55:00.000Z
        start_orbit_number: 102808
        end_orbit_number: 103246
        map_projection: Equidistant Cylindrical
        latitude_units: degrees_north
        longitude_units: degrees_east
        northernmost_latitude: 90
        southernmost_latitude: -90
        westernmost_longitude: -180
        easternmost_longitude: 180
        geospatial_lat_max: 90
        geospatial_lat_min: -90
        geospatial_lon_max: 180
        geospatial_lon_min: -180
        latitude_step: 0.0416666679084301
        longitude_step: 0.0416666679084301
        sw_point_latitude: -89.9791641235352
        sw_point_longitude: -179.97917175293
        spatialResolution: 4.64 km
        geospatial_lon_resolution: 0.0416666679084301
        geospatial_lat_resolution: 0.0416666679084301
        geospatial_lat_units: degrees_north
        geospatial_lon_units: degrees_east
        number_of_lines: 4320
        number_of_columns: 8640
        measure: Mean
        suggested_image_scaling_minimum: -2
        suggested_image_scaling_maximum: 45
        suggested_image_scaling_type: LINEAR
        suggested_image_scaling_applied: No
        _lastModified: 2021-12-03T08:21:22.000Z
        Conventions: CF-1.6 ACDD-1.3
        institution: NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group
        standard_name_vocabulary: CF Standard Name Table v36
        naming_authority: gov.nasa.gsfc.sci.oceandata
        id: AQUA_MODIS.20210901_20210930.L3b.MO.SST.nc/L3/AQUA_MODIS.20210901_20210930.L3b.MO.SST.nc
        license: https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/
        creator_name: NASA/GSFC/OBPG
        publisher_name: NASA/GSFC/OBPG
        creator_email: data@oceancolor.gsfc.nasa.gov
        publisher_email: data@oceancolor.gsfc.nasa.gov
        creator_url: https://oceandata.sci.gsfc.nasa.gov
        publisher_url: https://oceandata.sci.gsfc.nasa.gov
        processing_level: L3 Mapped
        cdm_data_type: grid
        keywords: Earth Science > Oceans > Ocean Optics > Sea Surface Temperature
        keywords_vocabulary: NASA Global Change Master Directory (GCMD) Science Keywords
        data_bins: 20227868
        data_minimum: -1.80000007152557
        data_maximum: 40.0000038146973
> SST_brick <- brick(list[59], "sst")
Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'brick': object of type 'builtin' is not subsettable
> ##Print the results for SST
> print(SST)
File /Users/kirstymedcalf/Documents/DMAD/Publication/DMAD_Maps_Analysis_Publication/Montenegro_Final_Analysis_Folders/Ocean_ColorSST_2016_2021/AQUA_MODIS.20210901_20210930.L3m.MO.SST.sst.4km.nc (NC_FORMAT_NETCDF4):

     3 variables (excluding dimension variables):
        short sst[lon,lat]   (Chunking: [87,44])  (Compression: shuffle,level 4)
            long_name: Sea Surface Temperature
            scale_factor: 0.00499999988824129
            add_offset: 0
            units: degree_C
            standard_name: sea_surface_temperature
            _FillValue: -32767
            valid_min: -1000
            valid_max: 10000
            display_scale: linear
            display_min: -2
            display_max: 45
        unsigned byte qual_sst[lon,lat]   (Chunking: [87,44])  (Compression: shuffle,level 4)
            long_name: Quality Levels, Sea Surface Temperature
            _FillValue: 255
            valid_min: 0
            valid_max: 5
        unsigned byte palette[eightbitcolor,rgb]   (Contiguous storage)  

     4 dimensions:
        lat  Size:4320 
            long_name: Latitude
            units: degrees_north
            standard_name: latitude
            _FillValue: -999
            valid_min: -90
            valid_max: 90
        lon  Size:8640 
            long_name: Longitude
            units: degrees_east
            standard_name: longitude
            _FillValue: -999
            valid_min: -180
            valid_max: 180
        rgb  Size:3 (no dimvar)
        eightbitcolor  Size:256 (no dimvar)
[1] ">>>> WARNING <<<  attribute data_bins is an 8-byte value, but R"
[1] "does not support this data type. I am returning a double precision"
[1] "floating point, but you must be aware that this could lose precision!"

    61 global attributes:
        product_name: AQUA_MODIS.20210901_20210930.L3m.MO.SST.sst.4km.nc
        instrument: MODIS
        title: MODISA Level-3 Standard Mapped Image
        project: Ocean Biology Processing Group (NASA/GSFC/OBPG)
        platform: Aqua
        temporal_range: month
        processing_version: R2019.0
        date_created: 2021-12-03T08:21:22.000Z
        history: l3mapgen par=AQUA_MODIS.20210901_20210930.L3m.MO.SST.sst.4km.nc.param 
        l2_flag_names: LAND,HISOLZEN
        time_coverage_start: 2021-09-01T00:45:00.000Z
        time_coverage_end: 2021-10-01T02:55:00.000Z
        start_orbit_number: 102808
        end_orbit_number: 103246
        map_projection: Equidistant Cylindrical
        latitude_units: degrees_north
        longitude_units: degrees_east
        northernmost_latitude: 90
        southernmost_latitude: -90
        westernmost_longitude: -180
        easternmost_longitude: 180
        geospatial_lat_max: 90
        geospatial_lat_min: -90
        geospatial_lon_max: 180
        geospatial_lon_min: -180
        latitude_step: 0.0416666679084301
        longitude_step: 0.0416666679084301
        sw_point_latitude: -89.9791641235352
        sw_point_longitude: -179.97917175293
        spatialResolution: 4.64 km
        geospatial_lon_resolution: 0.0416666679084301
        geospatial_lat_resolution: 0.0416666679084301
        geospatial_lat_units: degrees_north
        geospatial_lon_units: degrees_east
        number_of_lines: 4320
        number_of_columns: 8640
        measure: Mean
        suggested_image_scaling_minimum: -2
        suggested_image_scaling_maximum: 45
        suggested_image_scaling_type: LINEAR
        suggested_image_scaling_applied: No
        _lastModified: 2021-12-03T08:21:22.000Z
        Conventions: CF-1.6 ACDD-1.3
        institution: NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group
        standard_name_vocabulary: CF Standard Name Table v36
        naming_authority: gov.nasa.gsfc.sci.oceandata
        id: AQUA_MODIS.20210901_20210930.L3b.MO.SST.nc/L3/AQUA_MODIS.20210901_20210930.L3b.MO.SST.nc
        license: https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/
        creator_name: NASA/GSFC/OBPG
        publisher_name: NASA/GSFC/OBPG
        creator_email: data@oceancolor.gsfc.nasa.gov
        publisher_email: data@oceancolor.gsfc.nasa.gov
        creator_url: https://oceandata.sci.gsfc.nasa.gov
        publisher_url: https://oceandata.sci.gsfc.nasa.gov
        processing_level: L3 Mapped
        cdm_data_type: grid
        keywords: Earth Science > Oceans > Ocean Optics > Sea Surface Temperature
        keywords_vocabulary: NASA Global Change Master Directory (GCMD) Science Keywords
        data_bins: 20227868
        data_minimum: -1.80000007152557
        data_maximum: 40.0000038146973

##Extract SST, longitude, and latitude values from each Aqua Modis file to create 59 lists for the variables of interest. However, I think I am only extracting the values for a single Aqua Modis file [1]

SST_filenames <- ncvar_get(SST, "sst")
lon_filenames <- ncvar_get(SST, "lon")
lat_filenames <- ncvar_get(SST, "lat")

I thought I'd try the brick() function but unsuccessfully because I still think I'm extracting values from a single Aqua Modis file [1]. .

SST_brick <- brick(filenames[59], "sst")
lon_brick <- brick(filenames[59], "lon")
lat_brick <- brick(filenames[59], "lat")

Results for SST_brick <- brick(filenames[59], "sst"):

   lass      : RasterBrick 
    dimensions : 4320, 8640, 37324800, 1  (nrow, ncol, ncell, nlayers)
    resolution : 0.04166667, 0.04166667  (x, y)
    extent     : -180, 180, -90.00001, 90  (xmin, xmax, ymin, ymax)
    crs        : +proj=longlat +datum=WGS84 +no_defs 
    source     : AQUA_MODIS.20210901_20210930.L3m.MO.SST.sst.4km.nc 
    names      : layer 
    varname    : sst

If I use the object SST, I think I am still opening the values for a single Aqua Modis netCDF file [1]

SST_filenames <- ncvar_get(SST, "sst")
lon_filenames <- ncvar_get(SST, "lon")
lat_filenames <- ncvar_get(SST, "lat")
Alice Hobbs
  • 1,021
  • 1
  • 15
  • 31
  • If you have a spatial dataset, I wonder whether you could read in your dataset using the terra package, which can usually process spatial .nc files. You could read all files at once with something like files_path = list.files(); terra::rast(files_path), which will give you an object with as many layers as you have files, and then rasterize and do analysis from there. – tlhenvironment Apr 08 '22 at 08:44
  • I used list.files(), which has given me an object with many files and initially, I rasterized them using ncin_SST <- raster::stack(filenames, varname = "sst"). However, this code only pulls out the sea surface temperature values and not the units of longitude and latitude. Do you have any suggestions? I need the long & lat units for each raster file in the list to position the sea surface temperature values. Many thanks if you can help – Alice Hobbs Apr 08 '22 at 10:34
  • I definitely agree that the netCDF files need to be converted into an object with multiple raster files but I can't figure out how to do this and extract the SST, long, and lat variables beforehand. – Alice Hobbs Apr 08 '22 at 10:37

1 Answers1

2

You seem to misunderstand what a RasterStack or SpatRaster object is. You printed SST_brick and it clearly shows that it knows about the longitude and latitudes. You do not need to anything else.

With

library(terra)
filenames = list.files('~/Documents/Ocean_ColorSST_2016_2021',pattern='*.nc',full.names=TRUE)

And assuming that dolphins is a matrix or data.frame with variables "longitude" and "latitude", you could do this for one file like this:

SST <- rast(filenames[1], "sst")
e <- extract(SST, dolphins[, c("longitude", "latitude")])

But you can probably do this in one swoop like this

SSTs <- rast(filenames, "sst")
e <- extract(SSTs, dolphins[, c("longitude", "latitude")])
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Hello Robert. Thank you so much for answering this question. I have been so confused and I have been trying to figure out how to do this for almost 5 days. My shapefile is called 'points_spdf'. The object SSTs works great, and I am so happy. I assume 'dolphins' is my shape file. I ran e <- extract(SSTs, points_spdf[, c("longitude", "latitude")]) – Alice Hobbs Apr 08 '22 at 19:54
  • I got this error message: Error in h(simpleError(msg, call)) : error in evaluating the argument 'y' in selecting a method for function 'extract': argument "i" is missing, with no default In addition: Warning message: In points_spdf[, c("longitude", "latitude")] : j index ignored – Alice Hobbs Apr 08 '22 at 19:54
  • Please ask a separate question about extract if you cannot figure that out. It seems that you are using a `SpatialPoints` object. That would not match a SpatRaster but we do not know what exactly you are doing; making it difficult to help. – Robert Hijmans Apr 08 '22 at 20:30
  • Thanks for leading me down the right path. Part of my answer did not copy over. This is how I produced my shape file https://stackoverflow.com/questions/71763547/converting-the-coordinate-reference-system-crs-to-create-a-spatial-dataframe-u?noredirect=1#comment126858271_71763547 – Alice Hobbs Apr 08 '22 at 20:39
  • Hello Robert. I have tried to figure out how to extract the values, and I managed to do this successfully with another shapefile with the object SSTs and I extracted all the values per point. However, when I follow the same sequence with the data of interest, my output from the object using extract() are NAs. Would you be able to give me advice? https://stackoverflow.com/questions/71820236/extracting-raster-data-at-points-using-the-extract-function-how-to-prevent-th?noredirect=1#comment126918184_71820236. Many thanks in advance if this is possible – Alice Hobbs Apr 11 '22 at 04:20