2

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.

Improve this question
3
  • I'd argue that it's probably not too numerous/large to use gdalwarp, but if you have to do it in C++, try the warp api tutorial. You are mixing projected and geographic coordinate systems, and I suspect that is why some nodes aren't getting set, because they aren't there. Use the warp api and helper functions like GDALGetSuggestedWarpOutput.
    – user10353
    Dec 31, 2013 at 16:05
  • Kyle, thank you for your suggestion. My answer (below) is working. My solution uses GDALAutoCreateWarpedVRT(), which appears to use the warp API internally. I do wonder whether GDALAutoCreateWarpedVRT() is better or worse than using the warp API directly. Also, my solution is slower than I would like. I will try using gdal_translate to batch convert my source images into geographic projection.
    – Jason P
    Jan 8, 2014 at 19:11
  • AutoCreateWarpedVRT is a simplified implementation of the warp API. There is very little control over what it does, but if it works, it works. I think WarpOptions are passed through, and band mapping is supported, but little more.
    – user10353
    Jan 9, 2014 at 16:52

2 Answers 2

Reset to default
1

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
Improve this answer
7
  • 1
    Thank you for your help. I see from your explanation that there may be an implied scale factor of 111 kilometers/degree. But where does the GDAL API tell me the map projection's exact units and/or scale factor? I need programatic (API based) access to the map projection's actual units (output of the affine transformation). That .prj definition only mentions units of degrees, which doesn't make sense given the values of the affine transforms I am getting. Thanks!
    – Jason P
    Dec 20, 2013 at 16:35
  • In the prj, the last entry states that units for the projection are metres. If you use the sphere definition, the scale factor is exactly degree[Radians]*Radius.
    – AndreJ
    Dec 20, 2013 at 16:45
  • I'm sorry I missed the entry UNIT["Meter",1] in your comment. I added geographicRef.SetLinearUnits("Meter",1.0) and, which added UNIT["Meter",1] to my projection definition, and the affine transform remained unchanged. Now I need to determine the exact radius assumed by the equirectangular projection and verify the resulting lat/lons are correct. I did not see a sphere definition in gdal.org/ogr/classOGRSpatialReference.html So based on what you're saying plus the Equirectangular definition there should be a radius to convert meters to degrees or radians.
    – Jason P
    Dec 20, 2013 at 18:53
  • For the ellipsoid, the conversion would be the radius at the aequator, unless you specify the latitude of true scale lat_ts to another value than zero.
    – AndreJ
    Dec 21, 2013 at 6:18
  • Andre, Thanks again for your help. I've made some progress but still don't get the right results. I've edited my original post. Take a look and let me know.
    – Jason P
    Dec 21, 2013 at 21:15
1

My pixel measurements taken in GeoViewer were wrong because GeoViewer's pixel measurements are only correct at 100% zoom. That caused me to doubt the results I obtained yesterday.

So the answer is to use the "linear unit" scale factor of (6378137.0*3.14159265358979)/180.0: 

OGRSpatialReference geographicRef;
geographicRef.SetWellKnownGeogCS("WGS84");
geographicRef.SetProjection(SRS_PT_EQUIRECTANGULAR);
geographicRef.SetLinearUnits("My Degrees Lat,Lon",(6378137.0*3.14159265358979)/180.0);

This gives a geo transform matrix that produces latitudes and longitudes that differ by about 3e-7 degrees with respect to measurements taken in LizardTech GeoViewer and Google Maps. That's close enough. I also tested various UTM-projected GeoTiffs and got similar accuracy.

Question is now fully answered

Andre, thank you for your help!

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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