Fix CRS parser and proj4 cea projection support #2403

pomadchin · 2017-10-03T15:00:13Z

pomadchin · 2017-10-03T19:28:27Z

This PR fixes biggis-project/biggis-landuse#2 and cea projection support (before it produced black tile or even could not reproject anything). But this PR introduces a new question.

The rasters below:

(original in cea projection)

(reprojected into WebMercator via GDAL)

(reprojected into WebMercator via GeoTrellis)

The question is where is the mistake? In proj4j library (though i tried to check various of CylindricalEqualAreaProjection implementations in proj4, proj4js and jmaprojlib) or in GDAL (i don't think its a gdal issue though)?

For sure this bug (if it's a bug) can be addressed in a separate PR.

GDAL info:

Cea -> WebMercator (GeoTrellis):

$ gdalinfo -proj4 /tmp/cea_wm_gt.tiff 
Driver: GTiff/GeoTIFF
Files: /tmp/cea_wm_gt.tiff
Size is 352, 637
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
PROJ.4 string is:
'+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs'
Origin = (-13089969.861130239441991,5176329.847734357230365)
Pixel Size = (87.619885062502533,-87.619885062500600)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-13089969.861, 5176329.848) (117d35'21.12"W, 42d 6' 8.62"N)
Lower Left  (-13089969.861, 5120515.981) (117d35'21.12"W, 41d43'45.49"N)
Upper Right (-13059127.662, 5176329.848) (117d18'43.70"W, 42d 6' 8.62"N)
Lower Right (-13059127.662, 5120515.981) (117d18'43.70"W, 41d43'45.49"N)
Center      (-13074548.761, 5148422.914) (117d27' 2.41"W, 41d54'58.03"N)
Band 1 Block=352x12 Type=Byte, ColorInterp=Gray

Cea -> WebMercator (GDAL):

$ gdalinfo -proj4 ~/Downloads/cea_wm.tiff 
Driver: GTiff/GeoTIFF
Files: /Users/daunnc/Downloads/cea_wm.tiff
Size is 512, 517
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
PROJ.4 string is:
'+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs'
Origin = (-13095817.809482177719474,4021262.748792563565075)
Pixel Size = (72.328612721326948,-72.328612721326948)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-13095817.809, 4021262.749) (117d38'30.24"W, 33d56'37.80"N)
Lower Left  (-13095817.809, 3983868.856) (117d38'30.24"W, 33d39'52.95"N)
Upper Right (-13058785.560, 4021262.749) (117d18'32.64"W, 33d56'37.80"N)
Lower Right (-13058785.560, 3983868.856) (117d18'32.64"W, 33d39'52.95"N)
Center      (-13077301.685, 4002565.802) (117d28'31.44"W, 33d48'15.78"N)
Band 1 Block=512x16 Type=Byte, ColorInterp=Gray

jbouffard · 2017-10-03T20:43:40Z

The mistake appears to be with GeoTrellis. Using pyproj, this is the results that I get

In [21]: cea_box = box(-28493.167, 4255884.544, 2358.212, 4224973.143)

In [22]: project = partial(
    ...:     pyproj.transform,
    ...:     pyproj.Proj('+proj=cea +lat_ts=33.75 +lon_0=-117.333333333333 +x_0=0.0 +y_0=0.0 +datum=NAD27 +units=m '),
    ...:     pyproj.Proj('+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=
    ...: @null +wktext +no_defs'))

In [23]: transform(project, cea_box).bounds
Out[23]: 
(-13095817.8097412,
 3983903.7592389523,
 -13058750.785876077,
 4021262.7490324154)

This matches what GDAL returns, but differs in what GeoTrellis calculates.

pomadchin · 2017-10-03T22:12:09Z

I'll dig into it a bit more tomorrow

pomadchin · 2017-10-04T08:56:34Z

Fixed.

(reprojected into WebMercator via GeoTrellis)

GDALINFO:

$ gdalinfo -proj4 /tmp/cea_wm_gt.tif 
Driver: GTiff/GeoTIFF
Files: /tmp/cea_wm_gt.tif
Size is 512, 517
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
PROJ.4 string is:
'+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs'
Origin = (-13095719.172012956812978,4021260.549522797577083)
Pixel Size = (72.328088445188769,-72.328088445188129)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-13095719.172, 4021260.550) (117d38'27.05"W, 33d56'37.74"N)
Lower Left  (-13095719.172, 3983866.928) (117d38'27.05"W, 33d39'52.90"N)
Upper Right (-13058687.191, 4021260.550) (117d18'29.46"W, 33d56'37.74"N)
Lower Right (-13058687.191, 3983866.928) (117d18'29.46"W, 33d39'52.90"N)
Center      (-13077203.181, 4002563.739) (117d28'28.25"W, 33d48'15.73"N)
Band 1 Block=512x15 Type=Byte, ColorInterp=Gray

pomadchin · 2017-10-04T09:52:08Z

Fixed proj4j tests in proj4-epsg.csv (tests set to passing for all cea projections):

EPSG:4326→EPSG:3410
EPSG:4326→EPSG:3975
EPSG:4326→EPSG:6933

pomadchin force-pushed the fix/crs-parser branch from 4fa4274 to 0903246 Compare October 3, 2017 15:02

Fix our CRS Parser

8bf5051

pomadchin force-pushed the fix/crs-parser branch from 0903246 to 8bf5051 Compare October 3, 2017 15:29

pomadchin mentioned this pull request Oct 3, 2017

ERROR org.apache.spark.executor.Executor ... scala.MatchError: Some() at geotrellis.raster.io.geotiff.reader.GeoTiffCSParser.getEllipsoidInfo biggis-project/biggis-landuse#2

Closed

pomadchin requested a review from echeipesh October 3, 2017 15:31

pomadchin force-pushed the fix/crs-parser branch from 3f4f5c0 to c301d74 Compare October 3, 2017 17:44

Add projection support

a22042e

pomadchin force-pushed the fix/crs-parser branch from c301d74 to a22042e Compare October 3, 2017 17:45

More tests, fixed files

d341ed1

pomadchin force-pushed the fix/crs-parser branch from 29746ea to d341ed1 Compare October 3, 2017 18:17

echeipesh approved these changes Oct 3, 2017

View reviewed changes

pomadchin force-pushed the fix/crs-parser branch from 7c45c09 to bbdb6fa Compare October 4, 2017 08:59

Fix default trueScaleLatitude value for cea projection

b15961e

pomadchin force-pushed the fix/crs-parser branch from bbdb6fa to b15961e Compare October 4, 2017 09:04

pomadchin added 2 commits October 4, 2017 12:53

Set cea tests to passing

b4c3372

Added more lat/lon aliases

702d057

pomadchin force-pushed the fix/crs-parser branch from 6c34718 to 702d057 Compare October 4, 2017 10:10

pomadchin changed the title ~~Fix CRS Parser~~ Fix CRS parser and proj4 cea projection support Oct 4, 2017

pomadchin merged commit 6f1e0fb into locationtech:master Oct 4, 2017

pomadchin mentioned this pull request Dec 20, 2017

Problem with Austria_Lambert CRS #2528

Closed

pomadchin mentioned this pull request Mar 13, 2018

Backport Cylindrical Equal Area Projection fix locationtech/proj4j#8

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix CRS parser and proj4 cea projection support #2403

Fix CRS parser and proj4 cea projection support #2403

pomadchin commented Oct 3, 2017

pomadchin commented Oct 3, 2017 •

edited

Loading

jbouffard commented Oct 3, 2017

pomadchin commented Oct 3, 2017

pomadchin commented Oct 4, 2017 •

edited

Loading

pomadchin commented Oct 4, 2017 •

edited

Loading

Fix CRS parser and proj4 cea projection support #2403

Fix CRS parser and proj4 cea projection support #2403

Conversation

pomadchin commented Oct 3, 2017

pomadchin commented Oct 3, 2017 • edited Loading

jbouffard commented Oct 3, 2017

pomadchin commented Oct 3, 2017

pomadchin commented Oct 4, 2017 • edited Loading

pomadchin commented Oct 4, 2017 • edited Loading

pomadchin commented Oct 3, 2017 •

edited

Loading

pomadchin commented Oct 4, 2017 •

edited

Loading

pomadchin commented Oct 4, 2017 •

edited

Loading