GeoTiff writer does not currently support WebMercator projection with no EPSG code set

[wei-shujie]

https://github.com/locationtech/geotrellis/blob/b199f9474c59679a9ede074ea0774d973e424e81/raster/src/main/scala/geotrellis/raster/io/geotiff/writer/CoordinateSystemParser.scala#L159

Caused by: geotrellis.raster.io.geotiff.writer.GeoTiffWriterLimitationException: This GeoTiff writer does not currently support the projection +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs  without an EPSG code associated with the CRS. You'll need to use a CRS that has an EPSG code, or reproject before writing to GeoTIFF.
    at geotrellis.raster.io.geotiff.writer.CoordinateSystemParser.projProps$lzycompute(CoordinateSystemParser.scala:166)
    at geotrellis.raster.io.geotiff.writer.CoordinateSystemParser.projProps(CoordinateSystemParser.scala:158)
    at geotrellis.raster.io.geotiff.writer.CoordinateSystemParser.parse$lzycompute(CoordinateSystemParser.scala:129)
    at geotrellis.raster.io.geotiff.writer.CoordinateSystemParser.parse(CoordinateSystemParser.scala:99)
    at geotrellis.raster.io.geotiff.writer.CoordinateSystemParser$.parse(CoordinateSystemParser.scala:47)
    at geotrellis.raster.io.geotiff.writer.TiffTagFieldValue$.collect(TiffTagFieldValue.scala:149)
    at geotrellis.raster.io.geotiff.writer.GeoTiffWriter.$anonfun$appendCloudOptimized$1(GeoTiffWriter.scala:198)
    at scala.collection.immutable.List.map(List.scala:283)
    at geotrellis.raster.io.geotiff.writer.GeoTiffWriter.appendCloudOptimized(GeoTiffWriter.scala:197)
    at geotrellis.raster.io.geotiff.writer.GeoTiffWriter.write(GeoTiffWriter.scala:344)
    at geotrellis.raster.io.geotiff.writer.GeoTiffWriter$.write(GeoTiffWriter.scala:36)
    at geotrellis.raster.io.geotiff.GeoTiff.write(GeoTiff.scala:74)
    at geotrellis.raster.io.geotiff.GeoTiff.write$(GeoTiff.scala:73)
    at geotrellis.raster.io.geotiff.MultibandGeoTiff.write(MultibandGeoTiff.scala:26)
    at geotrellis.spark.store.file.cog.FileCOGLayerWriter.$anonfun$writeCOGLayer$2(FileCOGLayerWriter.scala:83)
    at geotrellis.spark.store.file.cog.FileCOGLayerWriter.$anonfun$writeCOGLayer$2$adapted(FileCOGLayerWriter.scala:78)

use FileCogLayerWriter to write and the crs is WebMercator will throw this exception，because the proj is merc，but CorrdinateSystemParser.scala Without this

[pomadchin]

Hi @iiuk3 , try to use geotrellis.proj4j.WebMercator projection to write TIFFs as a workaroud. It has an EPSG code properly set.

[pomadchin]

Also will leave a reference to params that have to be parsed https://github.com/OSGeo/gdal/blob/8d1cde33ca846c3a4071d5b7bf3035f2ad98e2f3/gdal/frmts/gtiff/gt_wkt_srs.cpp#L1701-L1733

[wei-shujie]

Hi @iiuk3 , try to use geotrellis.proj4j.WebMercator projection to write TIFFs as a workaroud. It has an EPSG code properly set.

But it still throws the exception this is my code

// a name for the layer to be ingested
    val layerName = "mod09a1_hdf"
    // the projection we'd like things to be tiled in
    val targetCRS = geotrellis.proj4.WebMercator
    // an interpolation method
    // the scheme for generating tile layouts
    val layoutScheme = ZoomedLayoutScheme(targetCRS, tileSize = 256)
    val paths = Seq("HDF4_EOS:EOS_GRID:\"Y:\\mod09a1\\2020.06.01\\MOD09A1.A2020153.h15v01.006.2020162062117.hdf\":MOD_Grid_500m_Surface_Reflectance:sur_refl_b01")
//    paths.map(uri => {
//      val metadata = GDALRasterSource(uri)
//      println(metadata)
//      metadata
//    })
    // Here, we parallelize a list of URIs and then turn them each into RasterSources
    // Note that reprojection/celltype manipulation is something that RasterSources allow us to do directly
    val sourceRDD: RDD[RasterSource] =
    sc.parallelize(paths, paths.length)
      .map(uri => GDALRasterSource(uri).reproject(targetCRS): RasterSource)
      .cache()
    // A RasterSummary is necessary for writing metadata; we can get one from an RDD
    val all = RasterSummary.collect(sourceRDD, _ => ())
    require(all.size == 1, "multiple CRSs detected")
    val summary = all.head

    // levelFor gives us the zoom and LayoutDefinition which most closely matches the computed Summary's Extent and CellSize
    val LayoutLevel(zoom, layout) = summary.levelFor(layoutScheme)

    // This is the actual in (spark, distributed) memory layer
    val contextRDD = RasterSourceRDD.tiledLayerRDD(sourceRDD, layout)

    // A reference to the attribute store for this layer
    val output = "D:/data/geotrellis"
    //val attributeStore = AttributeStore(output)

    // Build a layer writer
    val writer = FileCOGLayerWriter(output)

    // Actually write out the RDD constructed and transformed above
    writer.write(layerName,contextRDD,zoom,ZCurveKeyIndexMethod)

[pomadchin]

Hi @iiuk3 could you test your code with the following PR? https://github.com/locationtech/geotrellis/pull/3283

Thanks!

[pomadchin]

This issue was automatically closed, feel free to reopen it if the issue is still present in the current snapshot

[wei-shujie]

Hi @iiuk3 could you test your code with the following PR? #3283

Thanks!

I re-tested, this exception was resolved, but another exception occurred. The exception stack is as follows:

Exception in thread "main" java.util.NoSuchElementException: None.get
    at scala.None$.get(Option.scala:349)
    at scala.None$.get(Option.scala:347)
    at geotrellis.store.cog.vrt.VRT.toXML(VRT.scala:140)
    at geotrellis.store.cog.vrt.VRT.write(VRT.scala:150)
    at geotrellis.spark.store.file.cog.FileCOGLayerWriter.$anonfun$writeCOGLayer$1(FileCOGLayerWriter.scala:118)
    at geotrellis.spark.store.file.cog.FileCOGLayerWriter.$anonfun$writeCOGLayer$1$adapted(FileCOGLayerWriter.scala:69)
    at scala.collection.Iterator.foreach(Iterator.scala:944)
    at scala.collection.Iterator.foreach$(Iterator.scala:944)
    at scala.collection.AbstractIterator.foreach(Iterator.scala:1432)
    at scala.collection.IterableLike.foreach(IterableLike.scala:71)
    at scala.collection.IterableLike.foreach$(IterableLike.scala:70)
    at scala.collection.AbstractIterable.foreach(Iterable.scala:54)
    at geotrellis.spark.store.file.cog.FileCOGLayerWriter.writeCOGLayer(FileCOGLayerWriter.scala:69)
    at geotrellis.spark.store.cog.COGLayerWriter.write(COGLayerWriter.scala:98)
    at geotrellis.spark.store.cog.COGLayerWriter.write$(COGLayerWriter.scala:79)
    at geotrellis.spark.store.file.cog.FileCOGLayerWriter.write(FileCOGLayerWriter.scala:37)
    at geotrellis.spark.store.cog.COGLayerWriter.write(COGLayerWriter.scala:66)
    at geotrellis.spark.store.cog.COGLayerWriter.write$(COGLayerWriter.scala:58)
    at geotrellis.spark.store.file.cog.FileCOGLayerWriter.write(FileCOGLayerWriter.scala:37)

I tracked the exception and found that the error appeared here：

the proj4String is “+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs ” ，But the correct value should be “+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”.

[pomadchin]

Hey @iiuk3, I think it is more a PRO4J parser issue now;

val bad = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"

CRS.fromString(bad).toWKT // > None

Also this proj4 string comes from your TIFF metadata? How did you generate it.

Also using webmeractor explicitly should help:

val crs = WebMercator
val wktString = crs.toWKT() // > Some(PROJCS["WGS 84 / Pseudo-Mercator", GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]], PROJECTION["Popular Visualisation Pseudo Mercator", AUTHORITY["EPSG","1024"]], PARAMETER["semi_minor", 6378137.0], PARAMETER["latitude_of_origin", 0.0], PARAMETER["central_meridian", 0.0], PARAMETER["scale_factor", 1.0], PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], UNIT["m", 1.0], AXIS["Easting", EAST], AXIS["Northing", NORTH], AUTHORITY["EPSG","3857"]])

[pomadchin]

@iiuk3 I also created another issue to cover this problem: https://github.com/locationtech/proj4j/issues/65

[wei-shujie]

Also this proj4 string comes from your TIFF metadata? How did you generate it.

Also using webmeractor explicitly should help:

I generate it by GDALRasterSource(uri).reproject(WebMercator).crs.toProj4String. And, if I change the version of gdal-warp-bindings to 1.0.1 and gdal to 2.4, I can get a correct result

[pomadchin]

@iiuk3 It makes sense, probably the older version of GDAL had some bugs / features related to the CRS proj4 string representation.

I would recommend you to depend on the most up to date GeoTrellis (3.4.1), it already uses gdal-warp-bindings 1.0.1 and yep, using the most up to date GDAL 2.4.x versoin is also recommended.

Nice that it finally works :+1: