2

I am trying to reproject MODIS Level 3 chlorophyll data (for instance, this file) using gdal in R.

I used gdal_translate and gdalwarp functions. There was no crs information associated with A20121832012213.L3m_MO_CHL_chlor_a_4km.hdf file.

I used:

s_srs= "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +units=m +no_defs +R=6371007.181"

gdalwarp(source,dest,
s_srs="+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +units=m +no_defs +R=6371007.181",
t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",srcnodata=-32767,dstnodata=NA)

Although gdal_translate can translate hdf file to other formats, like img, this file cannot be reprojected by gdalwarp because there is an error:

GeoTransform values not available which is returned by function GetGeoTransform

PS:
Description of Level 3 data says:

This image is typically a scaled value (16bit), two-dimensional array of an Equidistant Cylindrical projection of the globe.

However MODIS Level 3 data is different from SeaWifs Level 3 data. There is no problem with reprojecting SeaWifs data, but not MODIS.

  • Which projection do you want? – user32309 Jul 18 '14 at 8:49
3

After all I found solution, but not sure if that is right and elegant.

I used proj4 string that SeaWifs level 3 have:

s_srs="+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

and overwrite georeferenced bounds by the values from SeaWifs hdf file

a_ullr =c(-20037508.343,10018754.171,20037508.343,-10018754.171) 

after this gdal_translate creates img with GeoTransform values

gdal_translate("A20121832012213.L3m_MO_CHL_chlor_a_4km.hdf", "A20121832012213.L3m_MO_CHL_chlor_a_4km.img", of="HFA",b=1,a_srs="+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs",a_ullr =c(-20037508.343,10018754.171,20037508.343,-10018754.171))

and gdalwarp can reproject this file

gdalwarp("A20121832012213.L3m_MO_CHL_chlor_a_4km.img","A20121832012213.L3m_MO_CHL_chlor_a_4km_reprojected.img",t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",srcnodata=-32767,dstnodata=NA)

Reprojection looks pretty accurate. If somebody can prove that this approach is acceptable?

  • You are basically assuming that your dataset uses the (deprecated) WGS 84 / Plate Carree projection (EPSG: 32662; epsg.io/32662), which is correct according to this document: oceancolor.gsfc.nasa.gov/DOCS/Ocean_Level-3_SMI_Products.pdf. Side remark: there apparently is some confusion about Plate Carree and World Equidistant Cylindrical. EPSG:4087 (epsg.io/4087) is the current code for WGS 84 / World Equidistant Cylindrical, but uses different projected bounds and would therefore result in a north-south shift when applied to the dataset in question. – Markus M. Jul 31 '15 at 1:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.