1

I am trying to download and plot the modis data from here http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.s.chl.a.sst.php

I regularly use .nc files from the ocean colour website but I can't seem to work with the hdf files the same way. It is possible to convert the hdf to a raster? 

gdalinfo("~vgpm.2008245.hdf")
gdal_translate(sds, dst_dataset = "~test.tiff")

sh: -c: line 0: unexpected EOF while looking for matching ``'
sh: -c: line 1: syntax error: unexpected end of file
NULL

rast <- raster("~test.tiff")

Error in .local(.Object, ...) :

Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", :
Cannot create a RasterLayer object from this file. (file does not exist)

I have tried the advice here Reading hdf files into R and converting them to geoTIFF rasters

and received this error: sh: -c: line 0: unexpected EOF while looking for matching ``' sh: -c: line 1: syntax error: unexpected end of file

I have also tried the suggestion here Converting HDF to georeferenced file (geotiff, shapefile) and the error is: Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file.

I just want to extract the data and plot the productivity data for a given location over time in a line plot alongside other data that I have...
thank you

liv
  • 93
  • 10
  • Duplicate of https://stackoverflow.com/questions/36772341/reading-hdf-files-into-r-and-converting-them-to-geotiff-rasters or https://stackoverflow.com/questions/35323495/how-to-make-a-rasterbrick-from-hdf5-files-r or https://stackoverflow.com/questions/48599875/converting-hdf-to-georeferenced-file-geotiff-shapefile or https://stackoverflow.com/questions/15974643/how-to-deal-with-hdf5-files-in-r – M-- May 13 '19 at 22:38

1 Answers1

0

On windows (and linux if the hdf drivers have been installed) you do not need to convert these files. You can use them directly 

library(raster)
x <- raster("vgpm.2008245.hdf")
x
#class      : RasterLayer 
#dimensions : 1080, 2160, 2332800  (nrow, ncol, ncell)
#resolution : 1, 1  (x, y)
#extent     : 0, 2160, 0, 1080  (xmin, xmax, ymin, ymax)
#crs        : NA 
#source     : vgpm.2008245.hdf 
#names      : vgpm.2008245

As you point out in the comments, the extent is not correct (and the CRS is not defined; and the NA flag is not set either). You can fix that like so. 

extent(x) <- extent(-180, 180, -90, 90)
crs(x) <- "+proj=longlat +datum=WGS84"
NAvalue(x) <- -9999

x
#class      : RasterLayer 
#dimensions : 1080, 2160, 2332800  (nrow, ncol, ncell)
#resolution : 0.1666667, 0.1666667  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#source     : vgpm.2008245.hdf 
#names      : vgpm.2008245

And if you wanted a tif you could then do

 r <- writeRaster(x, "vgpm.2008245.tif"

This did not work on windows in the past. And apparently it does not work on Mac out of the box --- depending how you install GDAL. I think you need to use homebrew to install gdal2 using "--with complete"

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • thank you @RobertHijmans. I set the work dir but library(raster) x <- raster("vgpm.2008245.hdf") x, gives an error Error in .local(.Object, ...) : Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. Is it because im using a Mac? I tried f <- "~/vgpm.2008245.hdf" file.exists(f) it says true so the file does exist at least! but I can't run writeRaster(x, "vgpm.2008245.tif" Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘writeRaster’ for signature ‘"character", "character" – liv May 14 '19 at 11:04
  • It is different (more difficult) on a Mac... – Robert Hijmans May 14 '19 at 15:34
  • Thank you. I downloaded windows parallel for Mac, downloaded r and ran the code and it worked. I noticed that the projection is different (extent 0, 2160, 0, 1080). I wasn't sure how to extract the data for a location so I converted the extent like what I had done before and extracted one location 40S, 180W. Unfortunately the results are absurd... am I transforming this incorrectly? library(raster) library(rgdal) x <- brick("vgpm.2008245.hdf") x extent(x) <- extent(0, 360, -90, 90) x[x < -9000] <- NA x <- rotate(x) x pnt = matrix(c(-180, -40), ncol = 2) pnt extract(x, pnt) – liv May 14 '19 at 18:25
  • I edited my response to address these issues. – Robert Hijmans May 15 '19 at 04:31