Converting HDF to georeferenced file (geotiff,

Question

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.

Answer 1

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:

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)