Opened 14 years ago

Closed 14 years ago

Last modified 13 years ago

#2706 closed defect (fixed)

Generalize, Clarify and Correct Equirectangular Projection

Reported by: warmerdam Owned by: warmerdam
Priority: high Milestone:
Component: OGR_SRS Version: svn-trunk
Severity: normal Keywords: equirectangular
Cc: dgrichard thare


It appears there is a mixup between different uses of the equirectangular projection with whether the latitude argument is a natural origin or a latitude of true scale (also often called standard parallel in this context).

In fact this projection can potentially have a origin and a latitude of true scale and they are generally distinct concepts.

It appears the old latitude_of_origin parameter is/was being mismapped to latitude of true scale in one or more of pds, isis2, isis3, HFA, usgs, and pci coordinate system handlers.

Didier has already extended things via SetEquirectangular2() to support both parameters, but there is some need to fix old drivers and converters that were wrong.

This ticket is derived from #2478, to capture the outstanding and difficult aspect.

Also, we should move to the name standard_parallel_1 instead of pseudo_standard_parallel_1 for equirectangular for similarity with ESRI and to avoid confusion.

Attachments (1)

fix_jp2.cpp (9.5 KB ) - added by warmerdam 13 years ago.
GeoJP2 CenterLat? to StdParallel?1 fixer utility (C++ source).

Download all attachments as: .zip

Change History (13)

comment:1 by warmerdam, 14 years ago

Status: newassigned

I have modified the isis, pds, hdf4, and hfa drivers to treat their lat parameters as standard_parallel_1 based on a review (r15824).

I have modified the geotiff driver to support saving and reading standard_parallel_1 in geotiff format (r15825).

I have modified pci, and usgs translation functions to treat lat as stdparallel (r15826).

I have modified various OGRSpatialReference stuff to use standard_parallel_1 instead of pseudo_standard_parallel_1 (r15827).

I have modified the proj.4 translation to treat origin and stdparallel1 distinctly and map to/from lat_0 and lat_ts respectively.

comment:2 by warmerdam, 14 years ago

(that last change was r15828).

comment:3 by warmerdam, 14 years ago

I'm going to add some autotest suite stuff tonight before closing this.

Also, someone ought to review the panorama, and ilwis converters as I suspect but cannot confirm that their equirectangular handling is wrong. It is likely not critical to fix this before 1.6 though.

comment:4 by Even Rouault, 14 years ago

Ok, with the whole patch series, for the dataset from #2478, the corner coordinates are identical to what they were in GDAL 1.5, so the regression in that case has gone.

Maybe we should advertize somehow that people who used previously the OSRSetEquirectangular() method will still get the same WKT than before, but not the same PROJ.4 string. So, it is possible that people using that method outside GDAL tree were using it because it gave the correct PROJ.4 string and didn't notice the wrong WKT (like our internal drivers - isis, etc - did before). If they are interested by getting the same PROJ.4 strings as before, they should migrate their code like that :

OSRSetEquirectangular(hSRS, lat, long, false_easting, false_northing) in GDAL < 1.6.0
OSRSetEquirectangular2(hSRS, 0, long, lat, false_easting, false_northing) in GDAL >= 1.6.0

comment:5 by Even Rouault, 14 years ago

I've translated the above dataset to GeoTIFF and read it back with gdalinfo and it appears that libgeotiff should be updated with the following patch to retrieve now the StdParallel1 param (actually, it's basically Didier's patch in, but with index 2 instead of 3.):

Index: frmts/gtiff/libgeotiff/geo_normalize.c
--- frmts/gtiff/libgeotiff/geo_normalize.c	(révision 15830)
+++ frmts/gtiff/libgeotiff/geo_normalize.c	(copie de travail)
@@ -1661,12 +1661,23 @@
                           &dfNatOriginLat, 0, 1 ) == 0 )
             dfNatOriginLat = 0.0;
+        if( psDefn->CTProjection == CT_Equirectangular )
+        {
+            if( GTIFKeyGet(psGTIF, ProjStdParallel1GeoKey, 
+                        &dfStdParallel1, 0, 1 ) == 0 )
+                dfStdParallel1 = 0.0;
+        }
         /* notdef: should transform to decimal degrees at this point */
         psDefn->ProjParm[0] = dfNatOriginLat;
         psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
         psDefn->ProjParm[1] = dfNatOriginLong;
         psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
+        if( psDefn->CTProjection == CT_Equirectangular )
+        {
+            psDefn->ProjParm[2] = dfStdParallel1;
+            psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
+        }
         psDefn->ProjParm[5] = dfFalseEasting;
         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
         psDefn->ProjParm[6] = dfFalseNorthing;

This also raises the question if we can deal with the fact that all existing geotiff (at least produced by GDAL) are using incorrectly the NatOriginLat instead of StdParallel1.

For fun, I've loaded my newly created geotiff with stdParallel1=60 in GlobalMapper10 and it reports 0 as the latitude of true scale, so many viewers around are likely confused too. should probably be updated to reflect the clarified meaning of parameters ?

comment:6 by warmerdam, 14 years ago

Resolution: fixed
Status: assignedclosed


I have patched libgeotiff similarly to what was suggested and pulled it back into GDAL (r15835).

I have also updated various test scripts (r15831, r15832)

I have modified ESRI translation to ensure that equirectangular maps to equidistant cylindrical, not plate carree. I have also stripped out the latitude_of_origin parameter in morphToESRI() if it is zero as the PE is not expecting it.

I also reordered the parameters in equirectangular so that stdpar1 comes before false easting/northing (r15834).

Now I can confirm that isis results appear to match the 1.5 corner coordinates with gdalinfo, and the tests seem to pass.

Didier, Trent, I'm going to prepare an RC1. I'd appreciate it if you would both carefully review equirectangular handling.

PS. I have rewritten the page.

comment:7 by warmerdam, 14 years ago

Incorporate the gt_srs_wkt.cpp changes into 1.6-esri (r17185).

in reply to:  6 comment:8 by dgrichard, 13 years ago

Didier, Trent, I'm going to prepare an RC1. I'd appreciate it if you would both carefully review equirectangular handling.

Made tests on today's trunk : it seems ok for me. I will have a look at the code.

comment:9 by thare, 13 years ago

Priority: highesthigh
Version: unspecifiedsvn-trunk

I have finally been able to take some time to look into these changes. I agree with the PDS/ISIS2/ISIS3 changes. (summary: all these readers define their latitude parameter as CENTER_LATITUDE, CENTER_LATITUDE, CenterLatitude respectively and those get mapped correctly to SetEquirectangular2 as the dfPseudoStdParallel1 parameter and NOT dfCenterLat parameter. Even though this is confusing it works.)

But I still have Equirectangular geoJpeg2000 Mars (HiRISE) issues. The Geotiff library "gt_wkt_srs.cpp" has also been changed to "SetEquirectangular2" but for these Mars images (which uses GDAL to help create the geotiff label for the jpeg2000 image - ) the "latitude" paramter is getting mapped to the dfCenterLat instead of the dfPseudoStdParallel1 parameter (see issue below).

I sort of understand the need to implement dfCenterLat and dfPseudoStdParallel1 but I'm not sure how we can be backwards compatible with already existing images. I guess I don't know where to go to help fix this? I don't think I can ask the HiRISE team to regenerate all their images (TBs worth - maybe eventually). Maybe a switch in gdal - but then apps using gdal in the future will not know to throw that. Maybe finding something special in these images but the Lunar, LROC Team, plan to create the same type of geoJpeg2000 images.

Any other ideas?

Here you can see the lat value is off the -+90 range (which breaks everything): ...

PROJECTIONEquirectangular, PARAMETER["latitude_of_origin",-45], PARAMETER["central_meridian",180], PARAMETER["standard_parallel_1",0], PARAMETER["false_easting",0], PARAMETER["false_northing",0], UNIT["metre",1,


Origin = (-6243330.250000000000000,-2815019.750000000000000) Pixel Size = (0.500000000000000,-0.500000000000000) Corner Coordinates: Upper Left (-6243330.250,-2815019.750) ( 74d21'32.27"E, 92d37'54.79"S) ...

Here is the image I have been testing - sorry it is so big.

thanks, Trent

comment:10 by warmerdam, 13 years ago

Milestone: 1.6.0


I have reviewed the file in question and it does appear the geotiff parameters are just wrong (done accordingly to the previous incorrect understanding). I don't see any general way that backward compatability can be offered.

What I could do is write a small utility that could be used to update the jpeg2000 files, in place, with the correct parameters. Thus, if the team producing the images can update their software to produce them properly, and if this utility was run over previously generated images we would hopefully be in a fairly good situation.

Let me know if you would like me to do so. The utility would just alter a few bytes in the geotiff section of the .jp2 files and would not require decompressing and recompressing the jpeg2000 imagery.

comment:11 by thare, 13 years ago

This all sounds good and I am comfortable with your assessment. I am glad you were willing to take one more look.

Quote: "I don't see any general way that backward compatibility can be offered." fair enough - let's move forward.

I would hate to ask for a work-around/patcher script but I think it could be useful for others. I also think, given an example, the code might be useful for other geoJpeg2000 tweaking tasks too.

I will request that the HiRISE Team to change/correct their code to meet with these findings. Fortunately this is a small change. By using another "trick" (e.g. setting a spherical body in the HiRISE label) we skirted the Equirectangular ocentric/ographic latitude system differences (the two systems define the scaling radius differently for this projection).

thanks, Trent

by warmerdam, 13 years ago

Attachment: fix_jp2.cpp added

GeoJP2 CenterLat? to StdParallel?1 fixer utility (C++ source).

comment:12 by warmerdam, 13 years ago

I have added a new version of the fix_jp2.cpp code that supports big endian or little endian geojp2 sections.

Note: See TracTickets for help on using tickets.