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?
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.
// 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.