Fix CRS parser and proj4 cea projection support by pomadchin · Pull Request #2403 · locationtech/geotrellis

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

Navigation Menu

Repository navigation

Fix CRS parser and proj4 cea projection support - #2403

Fix CRS parser and proj4 cea projection support#2403
pomadchin merged 6 commits into
locationtech:masterfrom
pomadchin:fix/crs-parser

pomadchin commented Oct 3, 2017

Uh oh!

pomadchin commented Oct 3, 2017 •

edited

Loading

Uh oh!

jbouffard commented Oct 3, 2017

Uh oh!

pomadchin commented Oct 3, 2017

Uh oh!

pomadchin commented Oct 4, 2017 •

edited

Loading

Uh oh!

pomadchin commented Oct 4, 2017 •

edited

Loading

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Conversation

pomadchin commented Oct 3, 2017

Uh oh!

pomadchin commented Oct 3, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jbouffard commented Oct 3, 2017

Uh oh!

pomadchin commented Oct 3, 2017

Uh oh!

pomadchin commented Oct 4, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pomadchin commented Oct 4, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

pomadchin commented Oct 3, 2017 •

edited

Loading

pomadchin commented Oct 4, 2017 •

edited

Loading

pomadchin commented Oct 4, 2017 •

edited

Loading