diff --git a/CHANGELOG.md b/CHANGELOG.md index f7f1c89326..f2602c4ec6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - GeoTrellisRasterSource should return None on empty reads [#3240](https://github.com/locationtech/geotrellis/pull/3240) - ArrayTile equals method always returns true if first elements are NaN [#3242](https://github.com/locationtech/geotrellis/issues/3242) - Fixed resource issue with JpegDecompressor that was causing a "too many open files in the system" exception on many parallel reads of JPEG compressed GeoTiffs. [#3249](https://github.com/locationtech/geotrellis/pull/3249) +- Fix MosaicRasterSource, GDALRasterSource and GeoTiffResampleRasterSource behavior [#3252](https://github.com/locationtech/geotrellis/pull/3252) ## [3.3.0] - 2020-04-07 diff --git a/gdal/src/main/scala/geotrellis/raster/gdal/GDALWarpOptions.scala b/gdal/src/main/scala/geotrellis/raster/gdal/GDALWarpOptions.scala index 010fa39776..2ff651ff55 100644 --- a/gdal/src/main/scala/geotrellis/raster/gdal/GDALWarpOptions.scala +++ b/gdal/src/main/scala/geotrellis/raster/gdal/GDALWarpOptions.scala @@ -256,20 +256,38 @@ case class GDALWarpOptions( */ def reproject(rasterExtent: GridExtent[Long], sourceCRS: CRS, targetCRS: CRS, resampleTarget: ResampleTarget = DefaultTarget, resampleMethod: ResampleMethod = NearestNeighbor): GDALWarpOptions = { val reprojectOptions = ResampleTarget.toReprojectOptions(rasterExtent, resampleTarget, resampleMethod) - val re = rasterExtent.reproject(sourceCRS, targetCRS, reprojectOptions) + val reprojectedRasterExtent = rasterExtent.reproject(sourceCRS, targetCRS, reprojectOptions) - this.copy( - cellSize = re.cellSize.some, - targetCRS = targetCRS.some, - sourceCRS = sourceCRS.some, - resampleMethod = reprojectOptions.method.some - ) + resampleTarget match { + case TargetDimensions(cols, rows) => + this.copy( + te = reprojectedRasterExtent.extent.some, + cellSize = None, + targetCRS = targetCRS.some, + sourceCRS = this.sourceCRS orElse sourceCRS.some, + resampleMethod = reprojectOptions.method.some, + dimensions = (cols.toInt, rows.toInt).some + ) + case _ => + val re = { + val targetRasterExtent = resampleTarget(reprojectedRasterExtent).toRasterExtent + if(this.alignTargetPixels) targetRasterExtent.alignTargetPixels else targetRasterExtent + } + + this.copy( + cellSize = re.cellSize.some, + te = re.extent.some, + targetCRS = targetCRS.some, + sourceCRS = this.sourceCRS orElse sourceCRS.some, + resampleMethod = reprojectOptions.method.some + ) + } } /** Adjust GDAL options to represents resampling with following parameters . * This call matches semantics and arguments of {@see RasterSource#resample} */ - def resample(gridExtent: => GridExtent[Long], resampleTarget: ResampleTarget): GDALWarpOptions = { + def resample(gridExtent: => GridExtent[Long], resampleTarget: ResampleTarget): GDALWarpOptions = resampleTarget match { case TargetDimensions(cols, rows) => this.copy(te = gridExtent.extent.some, cellSize = None, dimensions = (cols.toInt, rows.toInt).some) @@ -282,7 +300,7 @@ case class GDALWarpOptions( this.copy(te = re.extent.some, cellSize = re.cellSize.some) } - } + /** Adjust GDAL options to represents conversion to desired cell type. * This call matches semantics and arguments of {@see RasterSource#convert} */ diff --git a/gdal/src/test/resources/vlm/lc8-utm-1.tif b/gdal/src/test/resources/vlm/lc8-utm-1.tif new file mode 100644 index 0000000000..e73c65414c Binary files /dev/null and b/gdal/src/test/resources/vlm/lc8-utm-1.tif differ diff --git a/gdal/src/test/resources/vlm/lc8-utm-2.tif b/gdal/src/test/resources/vlm/lc8-utm-2.tif new file mode 100644 index 0000000000..3fa7f9cfc1 Binary files /dev/null and b/gdal/src/test/resources/vlm/lc8-utm-2.tif differ diff --git a/gdal/src/test/resources/vlm/lc8-utm-3.tif b/gdal/src/test/resources/vlm/lc8-utm-3.tif new file mode 100644 index 0000000000..d072b3ff92 Binary files /dev/null and b/gdal/src/test/resources/vlm/lc8-utm-3.tif differ diff --git a/gdal/src/test/resources/vlm/lc8-utm-re-mosaic-expected.tif b/gdal/src/test/resources/vlm/lc8-utm-re-mosaic-expected.tif new file mode 100644 index 0000000000..38d8e6a751 Binary files /dev/null and b/gdal/src/test/resources/vlm/lc8-utm-re-mosaic-expected.tif differ diff --git a/gdal/src/test/scala/geotrellis/raster/gdal/GDALRasterSourceSpec.scala b/gdal/src/test/scala/geotrellis/raster/gdal/GDALRasterSourceSpec.scala index ae09e84684..e235aac645 100644 --- a/gdal/src/test/scala/geotrellis/raster/gdal/GDALRasterSourceSpec.scala +++ b/gdal/src/test/scala/geotrellis/raster/gdal/GDALRasterSourceSpec.scala @@ -16,15 +16,15 @@ package geotrellis.raster.gdal -import geotrellis.proj4.LatLng +import cats.data.NonEmptyList +import geotrellis.proj4.{CRS, LatLng, WebMercator} import geotrellis.raster.geotiff.GeoTiffRasterSource import geotrellis.raster._ -import geotrellis.raster.io.geotiff.GeoTiff +import geotrellis.raster.io.geotiff.{AutoHigherResolution, GeoTiff} import geotrellis.raster.io.geotiff.reader.GeoTiffReader import geotrellis.raster.resample._ import geotrellis.raster.testkit._ import geotrellis.vector.Extent - import org.scalatest._ class GDALRasterSourceSpec extends FunSpec with RasterMatchers with GivenWhenThen { @@ -114,6 +114,45 @@ class GDALRasterSourceSpec extends FunSpec with RasterMatchers with GivenWhenThe source.metadata.attributes.mapValues(_.toUpperCase) shouldBe tsource.metadata.attributes.mapValues(_.toUpperCase) } + it("should perform a chained reprojection") { + val rs = GDALRasterSource(Resource.path("vlm/lc8-utm-1.tif")) + val drs = rs.reproject(WebMercator).reproject(LatLng).asInstanceOf[GDALRasterSource] + val rrs = rs.reproject(LatLng).asInstanceOf[GDALRasterSource] + + drs.crs shouldBe rrs.crs + + val RasterExtent(Extent(dxmin, dymin, dxmax, dymax), dcw, dch, dc, dr) = drs.gridExtent.toRasterExtent() + val RasterExtent(Extent(rxmin, rymin, rxmax, rymax), rcw, rch, rc, rr) = rrs.gridExtent.toRasterExtent() + + // our reprojection rounding error + // to fix it we can round cellSize in the options reprojection by 6 numbers after the decimal point + dxmin shouldBe rxmin +- 1e-1 + dymin shouldBe rymin +- 1e-1 + dxmax shouldBe rxmax +- 1e-1 + dymax shouldBe rymax +- 1e-1 + + dcw shouldBe rcw +- 1e-6 + dch shouldBe rch +- 1e-6 + + dc shouldBe rc +- 1 + dr shouldBe rr +- 1 + + // manually ensure that the difference in two GDALWarpOptions + drs.options.copy(cellSize = None, te = None) shouldBe rrs.options.copy(cellSize = None, te = None) + val CellSize(docw, doch) = drs.options.cellSize.get + val CellSize(rocw, roch) = rrs.options.cellSize.get + docw shouldBe rocw +- 1e-6 + doch shouldBe roch +- 1e-6 + + val Extent(dtexmin, dteymin, dtexmax, dteymax) = drs.options.te.get + val Extent(rtexmin, rteymin, rtexmax, rteymax) = rrs.options.te.get + + dtexmin shouldBe rtexmin +- 1e-1 + dteymin shouldBe rteymin +- 1e-1 + dtexmax shouldBe rtexmax +- 1e-1 + dteymax shouldBe rteymax +- 1e-1 + } + describe("should derive the cellType consistently with GeoTiffRasterSource") { val raster = Raster( ArrayTile(Array( @@ -158,5 +197,39 @@ class GDALRasterSourceSpec extends FunSpec with RasterMatchers with GivenWhenThe } } } + + + it("reprojectToRegion should produce an expected raster") { + val expected = GeoTiffRasterSource(Resource.path("vlm/lc8-utm-re-mosaic-expected.tif")).read().get + + val uris = List( + Resource.path("vlm/lc8-utm-1.tif"), + Resource.path("vlm/lc8-utm-2.tif"), + Resource.path("vlm/lc8-utm-3.tif") + ) + + val targetGridExtent = GridExtent[Long]( + Extent(-89.18876450968908, 37.796833507913995, -73.9486554460876, 41.41061537114088), + CellSize(0.01272458402544677,0.01272129304140357) + ) + val e = Extent(-89.18876450968908, 37.796833507913995, -73.9486554460876, 41.41061537114088) + val targetCRS = CRS.fromEpsgCode(4326) + + val rasterSources = NonEmptyList.fromListUnsafe(uris.map(GDALRasterSource(_))) + val moisac = MosaicRasterSource.instance(rasterSources, rasterSources.head.crs) + + val actual = + moisac + .reprojectToRegion( + targetCRS, + targetGridExtent.toRasterExtent, + NearestNeighbor, + AutoHigherResolution + ) + .read(e) + .get + + assertEqual(actual, expected) + } } } diff --git a/gdal/src/test/scala/geotrellis/raster/gdal/GDALWarpOptionsSpec.scala b/gdal/src/test/scala/geotrellis/raster/gdal/GDALWarpOptionsSpec.scala index 996ba6fb14..117405065e 100644 --- a/gdal/src/test/scala/geotrellis/raster/gdal/GDALWarpOptionsSpec.scala +++ b/gdal/src/test/scala/geotrellis/raster/gdal/GDALWarpOptionsSpec.scala @@ -23,6 +23,7 @@ import geotrellis.raster.resample._ import geotrellis.vector.Extent import geotrellis.raster.testkit._ +import cats.syntax.option._ import org.gdal.gdal._ import scala.collection.JavaConverters._ @@ -39,17 +40,19 @@ class GDALWarpOptionsSpec extends FunSpec with RasterMatchers with GivenWhenThen val reprojectOptions: GDALWarpOptions = generateWarpOptions( - sourceCRS = Some(CRS.fromString("+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +datum=NAD83 +units=m +no_defs ")), - targetCRS = Some(WebMercator), - cellSize = Some(CellSize(10, 10)) + sourceCRS = lccProjection.some, + targetCRS = WebMercator.some, + cellSize = CellSize(10, 10).some, + te = Extent(-8769160.0, 4257700.0, -8750640.0, 4274460.0).some ) val resampleOptions: GDALWarpOptions = generateWarpOptions( et = None, - cellSize = Some(CellSize(22, 22)), - sourceCRS = None, - targetCRS = None + cellSize = CellSize(22, 22).some, + sourceCRS = lccProjection.some, + targetCRS = WebMercator.some, + te = Extent(-8769178.0, 4257682.0, -8750640.0, 4274468.0).some ) def rasterSourceFromUriOptions(uri: String, options: GDALWarpOptions): GDALRasterSource = GDALRasterSource(uri, options) @@ -60,7 +63,7 @@ class GDALWarpOptionsSpec extends FunSpec with RasterMatchers with GivenWhenThen .EMPTY .reproject( rasterExtent = GridExtent(Extent(630000.0, 215000.0, 645000.0, 228500.0), 10, 10), - CRS.fromString("+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +datum=NAD83 +units=m +no_defs "), + lccProjection, WebMercator, TargetCellSize(CellSize(10, 10)) ) @@ -73,13 +76,13 @@ class GDALWarpOptionsSpec extends FunSpec with RasterMatchers with GivenWhenThen .EMPTY .reproject( rasterExtent = GridExtent(Extent(630000.0, 215000.0, 645000.0, 228500.0), 10, 10), - CRS.fromString("+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +datum=NAD83 +units=m +no_defs "), + lccProjection, WebMercator, TargetCellSize(CellSize(10, 10)) ) .resample( - GridExtent(Extent(-8769160.0, 4257700.0, -8750630.0, 4274460.0), CellSize(10, 10)), - TargetRegion(GridExtent[Long](Extent(-8769160.0, 4257700.0, -8750630.0, 4274460.0), CellSize(22, 22))) + GridExtent(Extent(-8769160.0, 4257700.0, -8750640.0, 4274460.0), CellSize(10, 10)), + TargetRegion(GridExtent[Long](Extent(-8769160.0, 4257700.0, -8750640.0, 4274460.0), CellSize(22, 22))) ) rasterSourceFromUriOptions(uri, opts) } @@ -141,12 +144,12 @@ class GDALWarpOptionsSpec extends FunSpec with RasterMatchers with GivenWhenThen val rs = GDALRasterSource(filePath) .reproject( - targetCRS = WebMercator, + targetCRS = WebMercator, resampleTarget = TargetCellSize(CellSize(10, 10)), - strategy = OverviewStrategy.DEFAULT + strategy = OverviewStrategy.DEFAULT ) .resampleToRegion( - region = GridExtent(Extent(-8769160.0, 4257700.0, -8750630.0, 4274460.0), CellSize(22, 22)) + region = GridExtent(Extent(-8769160.0, 4257700.0, -8750640.0, 4274460.0), CellSize(22, 22)) ) optimizedRawResample.gridExtent shouldBe rs.gridExtent @@ -164,11 +167,14 @@ class GDALWarpOptionsSpec extends FunSpec with RasterMatchers with GivenWhenThen } object GDALWarpOptionsSpec { + val lccProjection = CRS.fromString("+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +datum=NAD83 +units=m +no_defs ") + def generateWarpOptions( - et: Option[Double] = Some(0.125), - cellSize: Option[CellSize] = Some(CellSize(19.1, 19.1)), - sourceCRS: Option[CRS] = Some(CRS.fromString("+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +datum=NAD83 +units=m +no_defs ")), - targetCRS: Option[CRS] = Some(WebMercator) + et: Option[Double] = 0.125.some, + cellSize: Option[CellSize] = CellSize(19.1, 19.1).some, + sourceCRS: Option[CRS] = lccProjection.some, + targetCRS: Option[CRS] = WebMercator.some, + te: Option[Extent] = None ): GDALWarpOptions = { GDALWarpOptions( Some("VRT"), @@ -179,9 +185,9 @@ object GDALWarpOptionsSpec { None, sourceCRS, targetCRS, + te, None, - None, - List("-9999.0"), + Nil, Nil, Some(OverviewStrategy.DEFAULT), Nil, false, None, false, false, false, None, Nil, None, None, false, false, diff --git a/raster/data/vlm/lc8-utm-1-re.tif b/raster/data/vlm/lc8-utm-1-re.tif new file mode 100644 index 0000000000..d34377750e Binary files /dev/null and b/raster/data/vlm/lc8-utm-1-re.tif differ diff --git a/raster/data/vlm/lc8-utm-1.tif b/raster/data/vlm/lc8-utm-1.tif new file mode 100644 index 0000000000..e73c65414c Binary files /dev/null and b/raster/data/vlm/lc8-utm-1.tif differ diff --git a/raster/src/main/scala/geotrellis/raster/MosaicRasterSource.scala b/raster/src/main/scala/geotrellis/raster/MosaicRasterSource.scala index b2b31837eb..3b2787e207 100644 --- a/raster/src/main/scala/geotrellis/raster/MosaicRasterSource.scala +++ b/raster/src/main/scala/geotrellis/raster/MosaicRasterSource.scala @@ -19,8 +19,7 @@ package geotrellis.raster import geotrellis.vector._ import geotrellis.raster._ import geotrellis.raster.resample._ -import geotrellis.raster.reproject.Reproject -import geotrellis.proj4.{CRS, WebMercator} +import geotrellis.proj4.CRS import geotrellis.raster.io.geotiff.OverviewStrategy import cats.Semigroup @@ -86,13 +85,8 @@ abstract class MosaicRasterSource extends RasterSource { resampleTarget: ResampleTarget = DefaultTarget, method: ResampleMethod = ResampleMethod.DEFAULT, strategy: OverviewStrategy = OverviewStrategy.DEFAULT - ): RasterSource = - MosaicRasterSource( - sources map { _.reproject(targetCRS, resampleTarget, method, strategy) }, - crs, - gridExtent.reproject(this.crs, targetCRS, Reproject.Options.DEFAULT.copy(method = method)), - name - ) + ): RasterSource = + MosaicRasterSource.instance(sources.map(_.reproject(targetCRS, resampleTarget, method, strategy)), targetCRS, name) def read(extent: Extent, bands: Seq[Int]): Option[Raster[MultibandTile]] = { val rasters = sources map { _.read(extent, bands) } @@ -108,11 +102,13 @@ abstract class MosaicRasterSource extends RasterSource { resampleTarget: ResampleTarget, method: ResampleMethod, strategy: OverviewStrategy - ): RasterSource = MosaicRasterSource( - sources map { _.resample(resampleTarget, method, strategy) }, crs, name) + ): RasterSource = + MosaicRasterSource.instance(sources.map(_.resample(resampleTarget, method, strategy)), crs, name) def convert(targetCellType: TargetCellType): RasterSource = - MosaicRasterSource(sources map { _.convert(targetCellType) }, crs, name) + MosaicRasterSource.instance(sources.map(_.convert(targetCellType)), crs, name) + + override def toString: String = s"MosaicRasterSource(${sources.toList}, $crs, $gridExtent, $name)" } object MosaicRasterSource { @@ -123,7 +119,7 @@ object MosaicRasterSource { def combine(l: Raster[MultibandTile], r: Raster[MultibandTile]) = { val targetRE = RasterExtent( l.rasterExtent.extent combine r.rasterExtent.extent, - List(l.rasterExtent.cellSize, r.rasterExtent.cellSize).minBy(_.resolution)) + List(l.rasterExtent.cellSize, r.rasterExtent.cellSize).maxBy(_.resolution)) val result = l.resample(targetRE) merge r.resample(targetRE) result } @@ -144,13 +140,43 @@ object MosaicRasterSource { } } + /** + * This instance method allows to create an instance of the MosaicRasterSource. + * It assumes that all raster sources are known, and the GridExtent is also known. + */ + def instance(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS, targetGridExtent: GridExtent[Long], sourceName: SourceName): MosaicRasterSource = + new MosaicRasterSource { + val sources: NonEmptyList[RasterSource] = sourcesList + val crs: CRS = targetCRS + def gridExtent: GridExtent[Long] = targetGridExtent + val name: SourceName = sourceName + } + + /** + * This instance method allows to create an instance of the MosaicRasterSource. + * It computes the MosaicRasterSource basing on the input sourcesList. + */ + def instance(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS, sourceName: SourceName): MosaicRasterSource = { + val combinedExtent = sourcesList.map(_.extent).toList.reduce(_ combine _) + val minCellSize = sourcesList.map(_.cellSize).toList.maxBy(_.resolution) + val combinedGridExtent = GridExtent[Long](combinedExtent, minCellSize) + instance(sourcesList, targetCRS, combinedGridExtent, sourceName) + } + + def instance(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS): MosaicRasterSource = + instance(sourcesList, targetCRS, EmptyName) + + def instance(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS, targetGridExtent: GridExtent[Long]): MosaicRasterSource = + instance(sourcesList, targetCRS, targetGridExtent, EmptyName) + + /** All apply methods reproject the input sourcesList to the targetGridExtent */ def apply(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS, targetGridExtent: GridExtent[Long]): MosaicRasterSource = apply(sourcesList, targetCRS, targetGridExtent, EmptyName) def apply(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS, targetGridExtent: GridExtent[Long], rasterSourceName: SourceName): MosaicRasterSource = { new MosaicRasterSource { val name = rasterSourceName - val sources = sourcesList map { _.reprojectToGrid(targetCRS, gridExtent) } + val sources = sourcesList.map(_.reprojectToGrid(targetCRS, gridExtent)) val crs = targetCRS def gridExtent: GridExtent[Long] = targetGridExtent @@ -166,26 +192,11 @@ object MosaicRasterSource { val sources = sourcesList map { _.reprojectToGrid(targetCRS, sourcesList.head.gridExtent) } val crs = targetCRS def gridExtent: GridExtent[Long] = { - val reprojectedExtents = - sourcesList map { source => source.gridExtent.reproject(source.crs, targetCRS) } - val minCellSize: CellSize = reprojectedExtents.toList map { rasterExtent => - CellSize(rasterExtent.cellwidth, rasterExtent.cellheight) - } minBy { _.resolution } - reprojectedExtents.toList.reduce( - (re1: GridExtent[Long], re2: GridExtent[Long]) => { - re1.withResolution(minCellSize) combine re2.withResolution(minCellSize) - } - ) + val reprojectedSources = sources.toList + val combinedExtent = reprojectedSources.map(_.extent).reduce(_ combine _) + val minCellSize = reprojectedSources.map(_.cellSize).maxBy(_.resolution) + GridExtent[Long](combinedExtent, minCellSize) } } } - - @SuppressWarnings(Array("TraversableHead", "TraversableTail")) - def unsafeFromList(sourcesList: List[RasterSource], targetCRS: CRS = WebMercator, targetGridExtent: Option[GridExtent[Long]], rasterSourceName: SourceName = EmptyName): MosaicRasterSource = - new MosaicRasterSource { - val name = rasterSourceName - val sources = NonEmptyList(sourcesList.head, sourcesList.tail) - val crs = targetCRS - def gridExtent: GridExtent[Long] = targetGridExtent getOrElse { sourcesList.head.gridExtent} - } } diff --git a/raster/src/main/scala/geotrellis/raster/RasterSource.scala b/raster/src/main/scala/geotrellis/raster/RasterSource.scala index 578263baf8..7219e9605b 100644 --- a/raster/src/main/scala/geotrellis/raster/RasterSource.scala +++ b/raster/src/main/scala/geotrellis/raster/RasterSource.scala @@ -73,7 +73,7 @@ abstract class RasterSource extends CellGrid[Long] with RasterMetadata { */ def reprojectToRegion(targetCRS: CRS, region: RasterExtent, method: ResampleMethod = ResampleMethod.DEFAULT, strategy: OverviewStrategy = OverviewStrategy.DEFAULT): RasterSource = if (targetCRS == this.crs && region == this.gridExtent) this - else if (targetCRS == this.crs) resampleToRegion(region.asInstanceOf[GridExtent[Long]], method, strategy) + else if (targetCRS == this.crs) resampleToRegion(region.toGridType[Long], method, strategy) else reprojection(targetCRS, TargetRegion(region.toGridType[Long]), method, strategy) def resample(resampleTarget: ResampleTarget, method: ResampleMethod, strategy: OverviewStrategy): RasterSource diff --git a/raster/src/main/scala/geotrellis/raster/geotiff/GeoTiffReprojectRasterSource.scala b/raster/src/main/scala/geotrellis/raster/geotiff/GeoTiffReprojectRasterSource.scala index 49c1da2b07..09adb50ee2 100644 --- a/raster/src/main/scala/geotrellis/raster/geotiff/GeoTiffReprojectRasterSource.scala +++ b/raster/src/main/scala/geotrellis/raster/geotiff/GeoTiffReprojectRasterSource.scala @@ -71,13 +71,7 @@ class GeoTiffReprojectRasterSource( Reproject.Options.DEFAULT.copy(method = resampleMethod, errorThreshold = errorThreshold) ) - resampleTarget match { - case targetRegion: TargetRegion => targetRegion.region.toGridType[Long] - case targetAlignment: TargetAlignment => targetAlignment(reprojectedRasterExtent) - case targetDimensions: TargetDimensions => targetDimensions(reprojectedRasterExtent) - case targetCellSize: TargetCellSize => targetCellSize(reprojectedRasterExtent) - case _ => reprojectedRasterExtent - } + resampleTarget(reprojectedRasterExtent) } lazy val resolutions: List[CellSize] = diff --git a/raster/src/main/scala/geotrellis/raster/geotiff/GeoTiffResampleRasterSource.scala b/raster/src/main/scala/geotrellis/raster/geotiff/GeoTiffResampleRasterSource.scala index ca1ccc70a4..fee7eee963 100644 --- a/raster/src/main/scala/geotrellis/raster/geotiff/GeoTiffResampleRasterSource.scala +++ b/raster/src/main/scala/geotrellis/raster/geotiff/GeoTiffResampleRasterSource.scala @@ -74,13 +74,7 @@ class GeoTiffResampleRasterSource( Reproject.Options.DEFAULT.copy(method = resampleMethod, errorThreshold = errorThreshold) ) - resampleTarget match { - case targetRegion: TargetRegion => targetRegion.region.toGridType[Long] - case targetAlignment: TargetAlignment => targetAlignment(reprojectedRasterExtent) - case targetDimensions: TargetDimensions => targetDimensions(reprojectedRasterExtent) - case targetCellSize: TargetCellSize => targetCellSize(reprojectedRasterExtent) - case _ => reprojectedRasterExtent - } + resampleTarget(reprojectedRasterExtent) } } @@ -124,10 +118,7 @@ class GeoTiffResampleRasterSource( geoTiffTile.crop(windows.keys.toSeq, bands.toArray).map { case (gb, tile) => val targetRasterExtent = windows(gb) - Raster( - tile = tile, - extent = targetRasterExtent.extent - ).resample(targetRasterExtent.cols, targetRasterExtent.rows, method) + Raster(tile, closestTiffOverview.rasterExtent.extentFor(gb, clamp = true)).resample(targetRasterExtent, method) } } diff --git a/raster/src/test/scala/geotrellis/raster/MosaicRasterSourceSpec.scala b/raster/src/test/scala/geotrellis/raster/MosaicRasterSourceSpec.scala index 31653a32da..4cc1e81928 100644 --- a/raster/src/test/scala/geotrellis/raster/MosaicRasterSourceSpec.scala +++ b/raster/src/test/scala/geotrellis/raster/MosaicRasterSourceSpec.scala @@ -17,12 +17,14 @@ package geotrellis.raster import geotrellis.raster.geotiff._ -import geotrellis.raster.io.geotiff.GeoTiffTestUtils -import geotrellis.proj4.{LatLng, WebMercator} +import geotrellis.raster.io.geotiff.{AutoHigherResolution, GeoTiffTestUtils} +import geotrellis.proj4.{CRS, LatLng, WebMercator} import geotrellis.vector.Extent - import geotrellis.raster.testkit.RasterMatchers import cats.data.NonEmptyList +import cats.syntax.apply._ +import cats.instances.option._ +import geotrellis.raster.resample.NearestNeighbor import org.scalatest._ @@ -49,11 +51,11 @@ class MosaicRasterSourceSpec extends FunSpec with RasterMatchers with GeoTiffTes it("should union extents of its sources") { mosaicRasterSource.gridExtent shouldBe ( gtRasterSource1.gridExtent combine gtRasterSource2.gridExtent - ) + ) } it("should union extents with reprojection") { - mosaicRasterSource.reproject(WebMercator).gridExtent shouldBe mosaicRasterSource.gridExtent.reproject(LatLng, WebMercator) + mosaicRasterSource.reproject(WebMercator).gridExtent shouldBe mosaicRasterSource.sources.map(_.gridExtent.reproject(LatLng, WebMercator)).toList.reduce(_ combine _) } it("the extent read should match the extent requested") { @@ -88,10 +90,10 @@ class MosaicRasterSourceSpec extends FunSpec with RasterMatchers with GeoTiffTes val expectation = Raster( MultibandTile( IntConstantNoDataArrayTile(Array(1, 2, 3, 4, 1, 2, - 5, 6, 7, 8, 5, 6, - 9, 10, 11, 12, 9, 10, - 13, 14, 15, 16, 13, 14), - 6, 4)), + 5, 6, 7, 8, 5, 6, + 9, 10, 11, 12, 9, 10, + 13, 14, 15, 16, 13, 14), + 6, 4)), extentRead ) val result = mosaicRasterSource.read(extentRead, Seq(0)).get @@ -102,11 +104,11 @@ class MosaicRasterSourceSpec extends FunSpec with RasterMatchers with GeoTiffTes val expectation = Raster( MultibandTile( IntConstantNoDataArrayTile(Array(1, 2, 3, 4, 1, 2, 3, 4, - 5, 6, 7, 8, 5, 6, 7, 8, - 9, 10, 11, 12, 9, 10, 11, 12, - 13, 14, 15, 16, 13, 14, 15, 16), - 8, 4)), - mosaicRasterSource.gridExtent.extent + 5, 6, 7, 8, 5, 6, 7, 8, + 9, 10, 11, 12, 9, 10, 11, 12, + 13, 14, 15, 16, 13, 14, 15, 16), + 8, 4)), + mosaicRasterSource.gridExtent.extent ) val bounds = GridBounds(mosaicRasterSource.dimensions) val result = mosaicRasterSource.read(bounds, Seq(0)).get @@ -114,4 +116,47 @@ class MosaicRasterSourceSpec extends FunSpec with RasterMatchers with GeoTiffTes result.extent shouldEqual expectation.extent } } + + describe("reprojection operations") { + it("reprojectToRegion should behave consistent with simple raster sources") { + val rs = GeoTiffRasterSource(baseGeoTiffPath("vlm/lc8-utm-1.tif")) + val mrs = MosaicRasterSource(NonEmptyList(rs, Nil), rs.crs) + + rs.crs shouldBe mrs.crs + rs.gridExtent shouldBe mrs.gridExtent + + val targetCRS = CRS.fromEpsgCode(4326) + val targetGridExtent = GridExtent[Long]( + Extent(-76.66721364182322, 37.825520359396386, -73.9485516555973, 39.96447687368617), + CellSize(0.00355308391078038,0.00354453974736105) + ).toRasterExtent() + + val extent = Extent(-76.66721364182322, 37.825520359396386, -73.9485516555973, 39.96447687368617) + + val reprojectedrs = rs.reprojectToRegion( + targetCRS, + targetGridExtent, + NearestNeighbor, + AutoHigherResolution + ) + + val reprojectedmrs = mrs.reprojectToRegion( + targetCRS, + targetGridExtent, + NearestNeighbor, + AutoHigherResolution + ) + + reprojectedrs.crs shouldBe reprojectedmrs.crs + reprojectedrs.gridExtent shouldBe reprojectedmrs.gridExtent + + val result = reprojectedrs.read(extent) + val resultm = reprojectedmrs.read(extent) + + result shouldNot be (None) + resultm shouldNot be (None) + + (result, resultm).mapN(assertEqual) + } + } } diff --git a/raster/src/test/scala/geotrellis/raster/geotiff/GeoTiffRasterSourceSpec.scala b/raster/src/test/scala/geotrellis/raster/geotiff/GeoTiffRasterSourceSpec.scala index b9cf98daf2..1d7d94cd61 100644 --- a/raster/src/test/scala/geotrellis/raster/geotiff/GeoTiffRasterSourceSpec.scala +++ b/raster/src/test/scala/geotrellis/raster/geotiff/GeoTiffRasterSourceSpec.scala @@ -17,12 +17,11 @@ package geotrellis.raster.geotiff import geotrellis.raster._ -import geotrellis.raster.io.geotiff.GeoTiffTestUtils +import geotrellis.raster.io.geotiff.{AutoHigherResolution, GeoTiffTestUtils} import geotrellis.raster.io.geotiff.reader.GeoTiffReader import geotrellis.raster.resample._ import geotrellis.raster.testkit._ import geotrellis.vector._ - import org.scalatest._ class GeoTiffRasterSourceSpec extends FunSpec with RasterMatchers with GivenWhenThen with GeoTiffTestUtils { @@ -83,4 +82,21 @@ class GeoTiffRasterSourceSpec extends FunSpec with RasterMatchers with GivenWhen assertRastersEqual(actual, expected) } } + + it("resampleToRegion should produce an expected raster") { + val uri = baseGeoTiffPath("vlm/lc8-utm-1.tif") + val urie = baseGeoTiffPath("vlm/lc8-utm-1-re.tif") + + val expected = GeoTiffRasterSource(urie).read().get + + val re = GridExtent[Long](Extent(222285.0, 4187685.0, 589815.0, 4584315.0), CellSize(657.4776386404293,658.8538205980067)) + val e = Extent(222285.0, 4187685.0, 589815.0, 4584315.0) + + val actual = + GeoTiffRasterSource(uri) + .resampleToRegion(re, NearestNeighbor, AutoHigherResolution) + .read(e).get + + assertEqual(actual, expected) + } }