1

I am dealing with 28 HDF4 files on ocean primary productivity (annual .tar files can be found here: http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.cbpm2.v.php) My goal is to do some calculations (I need to calculate concentrations per area and obtain the mean over several years, i.e. combine all files spatially) and then convert them to a georeferenced file I can work with in ArcGIS (preferably shapefile, or geotiff).

I have tried several ways to convert to ASCII or raster files and then add a projection using gdalUtils tools such as gdal_translate and get_subdatasets. However, as the HDF4 files are not named after the standards (unlike the MODIS files), the latter doesn't work and I cannot access the subsets.

Here's the code I used to convert to raster:

library(raster)
library(gdalUtils)

setwd("...path_to_files...")

gdalinfo("cbpm.2015060.hdf")
hdf_file <- "cbpm.2015060.hdf"

outfile="testout"
gdal_translate(hdf_file,outfile,sds=TRUE,verbose=TRUE)
file.rename(outfile,paste("CBPM_test",".tif",sep="")) 

rast <- raster("CBPM_test.tif")

wgs1984 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
projection(rast) <- wgs1984
#crs(rast) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 

plot(rast)

writeRaster(rast, file="CBPM_geo.tif", format='GTiff', overwrite=TRUE)

The resulting projection is completely off. I'd appreciate help how to do this (converting through any format that works), preferably as batch process.

Kristina
  • 263
  • 2
  • 5
  • 11

1 Answers1

6

You've not set the extent of your raster, so its assumed to be 1:ncols, 1:nrows, and that's not right for a lat-long data set...

The gdalinfo implies its meant to be a full sphere, so if I do:

 extent(rast)=c(xmn=-180, xmx=180, ymn=-90, ymx=90)
 plot(rast)
 writeRaster(rast, "output.tif")

I see a raster with full global lat-long extent and when I load the raster into QGIS it overlays nicely with an OpenStreetMap.

There doesn't seem to be enough metadata in the file to do a precise projection (what's the Earth radius and eccentricity?) so don't try and do anything small-scale with this data...

Here's how it looks:

enter image description here

Also you've jumped through a few unnecessary hoops to read this. You can read the HDF directly and set its extent and projection:

> r = raster("./cbpm.2017001.hdf")

What do we get:

> r
class       : RasterLayer 
dimensions  : 1080, 2160, 2332800  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 2160, 0, 1080  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/rowlings/Downloads/HDF/cbpm.2017001.hdf 
names       : cbpm.2017001

Set extent: 

> extent(r)=c(xmn=-180, xmx=180, ymn=-90, ymx=90)

And projection:

> projection(r)="+init=epsg:4326"

And land values to NA:

> r[r==-9999]=NA

Write it, plot it:

> writeRaster(r,"r.tif")
> plot(r)
Spacedman
  • 92,590
  • 12
  • 140
  • 224
  • 1
    thank you very much for your help! Setting the extent did the trick indeed. However, trying to tidy up my, admittedly very messy, code, `r = raster("./cbpm.2017001.hdf")` gives me the error `Error in .local(.Object, ...) : Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file.` – Kristina Feb 03 '18 at 18:05
  • You'll have to do it your way - maybe your gdal doesn't have HDF drivers, but then gdal_translate works... hmmm sure that file exists for you? Its a different one than you used in your example. – Spacedman Feb 03 '18 at 18:43
  • Sorry for posting an example with another file but yes, the file exists and I cannot really figure out where the problem is. Would love to clean up my code as you suggested but for now it should work like it is. Thank you very much! Do you have any experience in performing calculations with raster files? This is what I am aiming at: http://orca.science.oregonstate.edu/faq01.php – Kristina Feb 03 '18 at 18:57
  • @Spacedman It must have something to do with fact that **raster** relies on **rgdal**, which comes with its own compiled version of **GDAL**, whereas **gdal_utils** uses system calls to whatever version of **GDAL** is on one's system path. (I seem to recall that rgdal, in turn, relies on a GDAL compiled by Brian Ripley.) Here on Windows, using the CRAN-distributed binaries for **raster** and **rgdal**, I get the same error as the OP. Not sure if this is relevant, but running `rgdal::gdalDrivers()` tells me that the HDF5 and HDF5Image drivers have values of `FALSE` in both `create` and `copy`. – Josh O'Brien Feb 03 '18 at 19:50
  • @Spacedman OK, that last bit of course **not** relevant. More to the point, `rgdal::gdalDrivers()` does not list any HDF4 driver. (And FWIW, the version of **GDAL** used on Windows by the **stars** package also does not support HDF4 (as shown in [this table](https://github.com/rwinlib/gdal2#full-config)).) – Josh O'Brien Feb 03 '18 at 20:08
  • FTR, thanks to requests documented [here](https://github.com/rwinlib/gdal2/issues/7), and Jeroen Ooms' generous work, it looks like `rgdal` and `sf` Windows binaries will soon ship with support for HDF4 drivers. – Josh O'Brien Feb 05 '18 at 19:04