Mosaic statistics in a custom projection

[bmcandr]

Hi,

I would like to calculate statistics for a mosaic in a custom projection (e.g., an equal area projection like EPSG:6933). I suspect this is not possible without modifications similar to those suggested in #591?

[vincentsarago]

Hi @bmcandr

In theory the mosaic backend feature() method accept both dst-crs and shape-crs options

titiler-pgstac/titiler/pgstac/mosaic.py

Lines 344 to 345 in d29a324

	dst_crs: Optional[CRS] = None,
	shape_crs: CRS = WGS84_CRS,

but they are not surfaced at the API level

titiler-pgstac/titiler/pgstac/factory.py

Lines 874 to 914 in d29a324

	def geojson_statistics(
	geojson: Union[FeatureCollection, Feature] = Body(
	..., description="GeoJSON Feature or FeatureCollection."
	),
	searchid=Depends(self.path_dependency),
	layer_params=Depends(self.layer_dependency),
	dataset_params=Depends(self.dataset_dependency),
	pixel_selection: PixelSelectionMethod = Query(
	PixelSelectionMethod.first, description="Pixel selection method."
	),
	max_size: int = Query(1024, description="Maximum image size to read onto."),
	stats_params=Depends(self.stats_dependency),
	histogram_params=Depends(self.histogram_dependency),
	pgstac_params: PgSTACParams = Depends(),
	backend_params=Depends(self.backend_dependency),
	reader_params=Depends(self.reader_dependency),
	env=Depends(self.environment_dependency),
	):
	"""Get Statistics from a geojson feature or featureCollection."""
	fc = geojson
	if isinstance(fc, Feature):
	fc = FeatureCollection(type="FeatureCollection", features=[geojson])
	threads = int(os.getenv("MOSAIC_CONCURRENCY", MAX_THREADS))
	with rasterio.Env(**env):
	with self.reader(
	searchid,
	reader_options={**reader_params},
	**backend_params,
	) as src_dst:
	for feature in fc:
	data, _ = src_dst.feature(
	feature.dict(exclude_none=True),
	pixel_selection=pixel_selection.method(),
	threads=threads,
	max_size=max_size,
	**layer_params,
	**dataset_params,
	**pgstac_params,
	)

2 solutions:

• add a dst-crs option
• add a shape-crs option but then you'll need to provide a geojson in a custom projection

both are available in titiler.core https://github.com/developmentseed/titiler/blob/main/src/titiler/core/titiler/core/dependencies.py#L454-L479

you can find examples of how to extend MosaicTilerFactories in https://stac-utils.github.io/titiler-pgstac/advanced/geojson_crop/ (note: the doc is not up to date and might try to import old functions)

[bmcandr]

Ah, I missed that DstCRSParams was already defined in the core. I had already implemented solution # 1 by adding near-identical dependency to the geojson_crop route and it works reasonably well.

That said, in order to ensure the cropped mosaic maintains the native resolution of the data I found that I had to expose width and height parameters as well. I then compute those based on the resolution of the underlying data and include them in the request to the tiler. Does that sound right?

PS Great job on this and all of the other tools you've shared with the community! Thanks for everything.

[vincentsarago]

Ah, I missed that DstCRSParams was already defined in the core. I had already implemented solution # 1 by adding near-identical dependency to the geojson_crop route and it works reasonably well.

it was added last week 🙈

That said, in order to ensure the cropped mosaic maintains the native resolution of the data I found that I had to expose width and height parameters as well. I then compute those based on the resolution of the underlying data and include them in the request to the tiler. Does that sound right?

Yeah it's always better to provide the height/width if you know it. I've always wondered about adding a resolution option but it will make things a bit more complex.

maybe it's because of the default max_size https://github.com/stac-utils/titiler-pgstac/blob/main/titiler/pgstac/mosaic.py#L346 which will be ignored when providing height/width

[bmcandr]

Well that explains why I didn't see DstCRSParams!

This is very helpful information. I'm glad to know I'm on the right track.