[GIS] GDAL/proj4 equirectangular projection coordinate system units assumed by OGRSpatialReference and GDALDataset::GetGeoTransform()

coordinate systemgdalprojraster

How do I configure my OGRSpatialReference object to specify an equirectangular projection system http://en.wikipedia.org/wiki/Equirectangular_projection where the affine geo-transform matrix (provided by GDALDataset::GetGeoTransform()) produces angular (WGS84 geodedic latitude & longitude) coordinates?

I did not find an answer here: http://www.gdal.org/gdal_datamodel.html or here: How do I convert affine coordinates to lat/lng? or in the GDALDataset and OGRSpatialReference class documentation.

My understanding of the equirectangular projection is that it is an angular (latitude & longitude) coordinate system by nature, possibly with scale factors involved. So given the affine transform matrix provided by GDALDataset::GetGeoTransform() for an equirectangular projected dataset, what is the nature of the projection's coordinate system? Is it angle-based (latitude & longitude) or based on virtual linear easting/northing distances that must then be converted to angles?

Note I am using a VRT dataset in order to query a UTM projected image as a virtual equirectangular (geographic) projected image.

My code:

// open a dataset that uses UTM projection, WGS84 coordinate system
if ((rawDataset = (GDALDataset *)GDALOpen(filename, GA_ReadOnly))==NULL) 
(throw exception)

(Output from rawDataset->GetProjectionRef(), etc. all look good. Confirm dataset uses UTM projection.)

 // specify a equirectangular projection (also known as geographic or uprojected lat/lon) for the VRT
OGRSpatialReference geographicRef;
geographicRef.SetProjection(SRS_PT_EQUIRECTANGULAR);
// geographicRef.SetEquirectangular(0,0,0,0); // gives same results as above
geographicRef.SetWellKnownGeogCS("WGS84"); // standard WGS84 coordinate system
geographicRef.SetProjCS("Equirectangular / WGS84");
//geographicRef.SetAngularUnits(SRS_UA_DEGREE,0.0174532925199433); // has no effect
//geographicRef.SetAngularUnits(SRS_UA_RADIAN,1.0); // has no effect 
geographicRef.exportToWkt(&pszDstWkt);

(Print the projection/coordinate specification string (pszDstWkt). Confirm projection is Equalrectangular and coordinate system is WGS84.

// create the VRT
if ((vrtGeographicDataset = (GDALDataset *)GDALAutoCreateWarpedVRT(rawDataset,(const char *)NULL,pszDstWkt,GRA_Bilinear,10,NULL))==NULL)
(throw exception)

Now examine the original image->UTM and VRT_image->equirectangular geo transform matrices:

rawDataset->GetGeoTransform(adfGeoTransform) gives values of:

387080.000000 0.100000 0.000000
3595399.000000 0.000000 -0.100000

The image has 0.1 meter pixel size. So the output of this transform is probably UTM easting/northing units of meters. However, it would be nice to find an API call that tells me the projection coordinate system units are meters.

However, vrtGeographicDataset->GetGeoTransform(adfGeoTransform) gives values of:

-11822338.617165 0.108005 0.000000
3616869.349192 0.000000 -0.108005

Clearly the X and Y offsets are too high to imply angular units of degrees or radians. And because 0.108005 is so close to 0.10000 I suspect the units are still proportional to meters. Either I am doing something wrong or perhaps there is an additional scaling factor necessary to convert to angles (degrees or radians).

Update on 2013-12-21

I've learned several things…

First GDALAutoCreateWarpedVRT() appears to be ignoring the destination coordinate system defined in char *pszDstWKT. OGRSpatialReference::SetGeogCS() and SetWellKnownGeogCS() have no effect. The absence of SetGeogCS() or SetWellKnownGeogCS() has no effect. Only SetProjection() and SetLinearUnits() appear to have an effect on the VRT's geo transform matrix.

Secondly, the VRT's geo transform is closer to what it should be relative to the source file's geo transform. The VRT geo transform has the correct signs while the original file's transform has signs consistent with the file's UTM projection. So GDALAutoCreateWarpedVRT() appears to be changing the projection to Equirectangular or something close to it.

Thirdly, via the following code I can obtain a geo transform that provides latitudes and longitudes that are somewhat close to expected values: 

OGRSpatialReference geographicRef;
geographicRef.SetProjection(SRS_PT_EQUIRECTANGULAR);
geographicRef.SetLinearUnits("My Degrees Lat,Lon",6378137.0*3.14159265358979/180.0); // equatorial radius
geographicRef.exportToWkt(&pszDstWkt);
vrtGeographicDataset = GDALAutoCreateWarpedVRT(rawDataset,NULL,pszDstWkt,...
vrtGeographicDataset->GetGeoTransform(adfGeoTransform);

Note pszDstWkt contains the definition: 'PROJCS["unnamed",PROJECTION["Equirectangular"],UNIT["My Degrees Lat,Lon",111319.4907932735]]'
The source file's definition is: 'PROJCS["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.2572235630016,

The VRT's geo transform is now:

-106.20187473836      9.70225388411317e-07                         0
32.4908901704265                         0     -9.70225388411317e-07

At pixel (0,0), the longitude of -106.20187473836 is exactly correct for my source image, verified by both the Lizardtech GeoViewer and Google Maps. However, the (0,0) pixel's latitude should be 32.490109, not 32.490890. And both diagonal coefficients of 9.70225388411317e-07 are somewhat too large. So all other pixels besides the upper left corner are in error.

Every other configuration of OGRSpatialReference I've tried gives worse results. And I don't see a more appropriate projection than equirectangular.

I did try to force a spherical definition (instead of ellipsoidal) by setting a high inverse flattening in geographicRef.SetGeogCS(). But again GDALAutoCreateWarpedVRT appears to be ignoring the PROJCS section of the definition altogether. So I am not sure how to achieve this.

Please help!! My goal is to extract raster sub-images (ROIs) with an affine geo transform that converts pixel (column,row) coordinates to WGS84 geodedic (longitude,latitude) coordinates. And I need to do this with GDAL's C++ API. The source files are too large and too numerous to use gdal_translate or gdalwarp.

Best Answer

To visualize the equirectangular projection, I chose this projection string :

+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs

and reprojected Natural Earths world boundaries and grid into it using QGIS:

enter image description here

The extent of the layers is +/- 20037400 horizontally and half of it vertically. So one degree is about 111km in both directions, what you might call a scale factor to transform directly from angular degrees to projected metres.

The WKT definition as stored in the .prj file is:

PROJCS["Equidistant_Cylindrical",
  GEOGCS["GCS_WGS_1984",
    DATUM["D_unknown",
    SPHEROID["WGS84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]],
  PROJECTION["Equidistant_Cylindrical"],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

If you don't like the results you get with that projection, try a sphere instead of the elliposid. For example EPSG:3786:

+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371007 +b=6371007 +units=m +no_defs
Related Question