TileLayerRDD.reproject error when RDD extent crosses antimeridian

[CloudNiner]

I attempted to reproject a TileLayerRDD[SpatialKey] with an extent that crosses the antimeridian because it contains tiled raster data for Fiji, whose extent also crosses the antimeridian.

Stacktrace:

 User class threw exception: org.apache.spark.SparkException: Job aborted due to stage failure: Task 3 in stage 9909.0 failed 4 times, most recent failure: Lost task 3.3 in stage 9909.0 (TID 309185, ip-172-31-62-66.ec2.internal, executor 39): java.lang.IllegalArgumentException: Encountered NaN during a refinement step: (NaN / 11827.774292424321). Input Extent(179.89374929576, -12.674583081471548, 180.10708262824002, -12.461249749) is likely not in source projection.
at geotrellis.vector.reproject.Reproject$.refine$1(Reproject.scala:104)
at geotrellis.vector.reproject.Reproject$.reprojectExtentAsPolygon(Reproject.scala:115)
at geotrellis.vector.reproject.Implicits$ReprojectExtent.reprojectAsPolygon(Implicits.scala:61)
at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:46)
at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:91)
at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:94)
at geotrellis.spark.reproject.TileRDDReproject$$anonfun$matchReprojectRasterExtent$1$$anonfun$14$$anonfun$apply$4$$anonfun$apply$5.apply(TileRDDReproject.scala:368)
at geotrellis.spark.reproject.TileRDDReproject$$anonfun$matchReprojectRasterExtent$1$$anonfun$14$$anonfun$apply$4$$anonfun$apply$5.apply(TileRDDReproject.scala:366) 

[pomadchin]

@CloudNiner could be helpful if you could document what is get passed into extent.reprojectAsPolygon(transform, 0.001).extent, otherwise it it hardly possible to track down the issue. Usually it meant that smth was wrong in the (extent, crs) tuple - extent didn't match the CRS (for instant extent is in some UTM zone, but GeoTrellis reader reads it in as a LatLng coordinates and tries to perform an incorrect transformation).

[echeipesh]

scala> val re = RasterExtent( Extent(179.89374929576, -12.674583081471548, 180.10708262824002, -12.461249749), 256,256)
re: geotrellis.raster.RasterExtent = GridExtent(Extent(179.89374929576, -12.674583081471548, 180.10708262824002, -12.461249749),CellSize(8.3333333000013E-4,8.333333299669829E-4),256x256)

scala> re.reproject(LatLng, WebMercator)
java.lang.IllegalArgumentException: Encountered NaN during a refinement step: (NaN / 11827.774292424321). Input Extent(179.89374929576, -12.674583081471548, 180.10708262824002, -12.461249749) is likely not in source projection.
  at geotrellis.vector.reproject.Reproject$.refine$1(Reproject.scala:104)
  at geotrellis.vector.reproject.Reproject$.reprojectExtentAsPolygon(Reproject.scala:115)
  at geotrellis.vector.reproject.Implicits$ReprojectExtent.reprojectAsPolygon(Implicits.scala:61)
  at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:46)
  at geotrellis.raster.GridExtent$gridExtentMethods$$anonfun$reproject$2.apply(GridExtent.scala:428)
  at geotrellis.raster.GridExtent$gridExtentMethods$$anonfun$reproject$2.apply(GridExtent.scala:428)
  at scala.Option.getOrElse(Option.scala:121)
  at geotrellis.raster.GridExtent$gridExtentMethods.reproject(GridExtent.scala:428)
  at geotrellis.raster.GridExtent$gridExtentMethods.reproject(GridExtent.scala:431)
  ... 40 elided

scala> re.extent.toPolygon.toGeoJson
res2: String = {"type":"Polygon","coordinates":[[[179.89374929576,-12.674583081471548],[179.89374929576,-12.461249749],[180.10708262824002,-12.461249749],[180.10708262824002,-12.674583081471548],[179.89374929576,-12.674583081471548]]]}

[pomadchin]

@echeipesh but @CloudNiner had a different exception O: Yours looks much better (at least it is not an infinite recursion and definitely a bug).

[echeipesh]

So what is happening here is that one of the x points is over 180 degrees. Proj4j transformation clips any x greater than 180 degrees to 180:

// this is x2
scala> JTS.Point(180.10708262824002, -12.674583081471548).reproject(LatLng, WebMercator)
res3: org.locationtech.jts.geom.Point = POINT (20037508.342789244 -1422578.3072905561)

// this is midpoint
scala> JTS.Point(180.000415962, -12.674583081471548).reproject(LatLng, WebMercator)
res4: org.locationtech.jts.geom.Point = POINT (20037508.342789244 -1422578.3072905561)

// this is just exactly 180
scala> JTS.Point(180, -12.674583081471548).reproject(LatLng, WebMercator)
res7: org.locationtech.jts.geom.Point = POINT (20037508.342789244 -1422578.3072905561)

// stepping down  by .1 to see change in projected X
scala> JTS.Point(179.900415962, -12.674583081471548).reproject(LatLng, WebMercator)
res5: org.locationtech.jts.geom.Point = POINT (20026422.698387947 -1422578.3072905561)

// stepping down  by .2 to see change in projected X
scala> JTS.Point(179.800415962, -12.674583081471548).reproject(LatLng, WebMercator)
res6: org.locationtech.jts.geom.Point = POINT (20015290.749308616 -1422578.3072905561)

[pomadchin]

ahhh I see; right, @CloudNiner has the same problem with the incorrect input?

[echeipesh]

I'm not sure, the symptom is that midpoint and one of the endpoints that are pest meridian are going to look the same, so value of deflect will be result of division by 0, Double.NaN.

If you look carefully at the stack trace in the issue it happens on the very first invocation of refine function, so there is no recursion happening, its straight to error.

[pomadchin]

Good catch, nice. I think this is probably not a bug though... What can be a proper label for such issues?

[echeipesh]

Well, its certainly not helpful behavior. I'm still not clear about the root cause. Diving into proj4j code there is nothing that would explain the clipping in projection:

https://github.com/locationtech/proj4j/blob/5ba8168a5e0346b79e5f7d2fd2bc25a68fba0670/src/main/java/org/locationtech/proj4j/proj/MercatorProjection.java#L32-L41

[echeipesh]

Related: https://github.com/locationtech/geotrellis/issues/2913

[echeipesh]

There is probably a whole set of issues berried here. This walks down the line of looking at https://github.com/locationtech/spatial4j integration or perhaps another method from PROJ upstream I'm not aware of

There is application level solution to THIS problem.

Since we're re projecting to a ZXY layout in WebMercator where tiles never cross anti-meridian we can fix the problem on the input side: don't reproject tiles that intersect anti-merdian, cut them before and reproject two halves. (I hope the pixels don't cross this line?)

That will leave the input valid and we know that each output tile will be valid by definition.

[jdries]

I also see this issue when doing: GDALRasterSource(path).reproject(crs) Where crs is Web Mercator and path is a global geotiff in EPSG:4326 with Extent(-180.0, -90.0, 180.0, 90.0).

So I get this error:

java.lang.IllegalArgumentException: Encountered NaN during a refinement step: (NaN / NaN). Input Extent(-180.0, -90.0, 180.0, 90.0) is likely not in source projection.
    at geotrellis.vector.reproject.Reproject$.refine$1(Reproject.scala:104)
    at geotrellis.vector.reproject.Reproject$.reprojectExtentAsPolygon(Reproject.scala:115)
    at geotrellis.vector.reproject.Implicits$ReprojectExtent.reprojectAsPolygon(Implicits.scala:61)
    at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:46)
    at geotrellis.raster.GridExtent$gridExtentMethods$$anonfun$reproject$2.apply(GridExtent.scala:428)
    at geotrellis.raster.GridExtent$gridExtentMethods$$anonfun$reproject$2.apply(GridExtent.scala:428)
    at scala.Option.getOrElse(Option.scala:121)
    at geotrellis.raster.GridExtent$gridExtentMethods.reproject(GridExtent.scala:428)
    at geotrellis.raster.gdal.GDALWarpOptions.reproject(GDALWarpOptions.scala:259)
    at geotrellis.raster.gdal.GDALRasterSource.reprojection(GDALRasterSource.scala:114)
    at geotrellis.raster.RasterSource.reproject(RasterSource.scala:56)

The problem here is that +-90 is beyond the valid extent of web mercator. However, only -90 seems to be a problem, as shown by this snippet:

    val t = Transform(LatLng, WebMercator)
    println(t.apply(180,90))
    println(t.apply(180,-90))

prints: (2.0037508342789244E7,2.3810769326496765E8) (2.0037508342789244E7,-Infinity) It is the -Infinity that triggers the NaN later on.

The consequence of this is that a valid global geotiff triggers this error when trying to convert to a webmercator tilesource.

[lukestewart13]

I also encountered this issue while attempting to reproject CHELSA climate data geotiffs from 4326 to 3857 on a RDD[SpatialKey, Tile] with Metadata[TileLayerMetadata[SpatialKey]]. In particular, this step

sourceRDD.reproject(WebMercator, layoutDefinition, Reproject.Options(...))

yielded the error:

IllegalArgumentException: Encountered NaN during a refinement step: (NaN / NaN). Input Extent(-180.00013888885002, -90.00847222214989, 180.00819300444988, 83.99986041515001) is likely not in source projection

The gdalinfo for a typical file looks like this (note that the geotiff extends beyond -180 and -90):

Driver: GTiff/GeoTIFF
Files: input.tif
Size is 43200, 20880
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000138888850017,83.999860415150010)
Pixel Size = (0.008333333300000,-0.008333333300000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0001389,  83.9998604) (180d 0' 0.50"W, 83d59'59.50"N)
Lower Left  (-180.0001389, -90.0001389) (180d 0' 0.50"W, 90d 0' 0.50"S)
Upper Right ( 179.9998597,  83.9998604) (179d59'59.49"E, 83d59'59.50"N)
Lower Right ( 179.9998597, -90.0001389) (179d59'59.49"E, 90d 0' 0.50"S)
Center      (  -0.0001396,  -3.0001392) (  0d 0' 0.50"W,  3d 0' 0.50"S)
Band 1 Block=43200x1 Type=UInt16, ColorInterp=Gray

Attempts to window the sourced geotiff failed. For example,

val geometry: Geometry = Polygon((-180.0, -90.0), (-180.0, 90.0), (180.0, 90.0), (180.0, -90.0), (-180.0, -90.0))

val rdd = HadoopGeoTiffRDD[ProjectedExtent, ProjectedExtent, Tile](
  path = new Path(...),
  uriToKey = (_, key) => key,
  options = Options.DEFAULT,
  geometry = Some(geometry)
)

yielded the same error, and it wasn't until I shrank the bounds to an intolerable extend (e.g. -175, -85, ...) that the job would succeed. Likewise, using gdalwarp to preprocess the file with -te set to various bounds also failed until I'd shrunk them to an intolerable extent.

The workaround was to use gdalwarp to do the reprojection from 4326 to 3857.