locationtech / geotrellis

GeoTrellis is a geographic data processing engine for high performance applications.
http://geotrellis.io
Other
1.34k stars 360 forks source link

Reproject and pyramid results in 'stretched' pixels #3206

Open marygriffus opened 4 years ago

marygriffus commented 4 years ago

The pyramids generated by the following result in 'stretched' pixels, with a pixel size of (1.396611782340103,-35.311459771682358).

      val layoutScheme = CustomZoomedLayoutScheme(crs, tileSize) // This is to provide the WorldExtent manually for our CRS

      val LayoutLevel(_, layoutDefinition) = layoutScheme.levelForZoom(maxZoom)
      val (_, retiled) = inputRDD.reproject(crs, layoutDefinition, Options(targetCellSize = Some(layoutDefinition.cellSize)))

      Pyramid.levelStream(retiled, layoutScheme, maxZoom, minZoom, Bilinear)
        .foreach {
          case (z, levelLayer) => {
            val keyedPath = "%s/%d".format(fileBase, z)
            val keyDir = new File(keyedPath)
            keyDir.mkdir()

            levelLayer.foreach({
              case (key, tile) => {
                val extent = levelLayer.metadata.layout.mapTransform(key)
                val gTiff = GeoTiff(tile, extent, levelLayer.metadata.crs)
                val path = s"$keyedPath/${fileName}_${key.col}_${key.row}.tif"
                GeoTiffWriter.write(gTiff, path)
              }
            })
          }
        }

Looks pretty normal zoomed out:

Screen Shot 2020-03-23 at 11 31 02 AM

But zoom in and the distortion is more clear:

Screen Shot 2020-03-23 at 11 32 01 AM Screen Shot 2020-03-23 at 11 31 21 AM
pomadchin commented 4 years ago

Hi @marygriffus, thanks for reporting! I marked it as a bug for now, probably requires some investigation. Could you provide some input raster to reproduce that behavior?

marygriffus commented 4 years ago

Will do!

marygriffus commented 4 years ago

Ok, here's a simplified version of file read/initial reprojection:

    val floatingLayout = FloatingLayoutScheme(512)

    val options: HadoopGeoTiffRDD.Options = HadoopGeoTiffRDD.Options(maxTileSize = Some(512))
    val baseRDD: RDD[(ProjectedExtent, Tile)] = HadoopGeoTiffRDD.spatial(tileDirectory, options)
    val meta: TileLayerMetadata[SpatialKey] = CollectTileLayerMetadata.fromRDD(baseRDD, floatingLayout)._2
    val gridded: RDD[(SpatialKey, Tile)] = baseRDD.tileToLayout(meta, NearestNeighbor)
    val inputRDD: TileLayerRDD[SpatialKey] = new ContextRDD(gridded, meta)

And here is a tile to replicate the issue. Mannassas_10m.tif.zip