OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.96k stars 2.58k forks source link

gdalwarp odd antimeridian handling #11400

Open Plantain opened 1 week ago

Plantain commented 1 week ago

What is the bug?

I have a space-projected GeoTIFF that spans the antimeridian - ( http://static.skysight.io/raw.tiff )

I want to transform it to EPSG:3857, and clip only a square from the "middle" of the data by defining extents that cross the antimeridian. I am not sure yet if it matters if the resulting image wraps or not, whether it could be 0-360 or must be -180 to 180, the result will have overviews generated and be used for tile serving. I am trying to avoid approaches that might lose fidelity or add processing time as this is to be done in realtime.

The image looks something like this: Screenshot 2024-12-01 at 11 29 13 am And I want to extract the middle: Screenshot 2024-12-01 at 11 29 32 am

I tried several approaches, and I am not sure if I am bouncing off bugs, undefined behaviour, or my own incompetence.

Steps to reproduce the issue

Firstly, warping to gdalwarp -of PNG -t_srs EPSG:3857 raw.tiff out.png Works as expected Screenshot 2024-12-01 at 11 11 17 am

Attempting to crop using extents beyond -180 crashes, I guess to be expected: gdalwarp -of PNG -t_srs EPSG:3857 raw.tiff out.png -te_srs EPSG:4326 -190 -65 -90 70 ERROR 1: Invalid dfWestLongitudeDeg Creating output file that is -39965P x 28552L. ERROR 1: Attempt to create -39965x28552 dataset is illegal,sizes must be larger than zero.

Interesting, the same in EPSG:3857 coordinates does run! But the result is missing everything west of the antimeridian gdalwarp -of PNG -t_srs EPSG:3857 raw.tiff out.png -te -22094105.797 -11068715.659 -10018754.171 11068715.659 ERROR 1: Point outside of projection domain Creating output file that is 23198P x 42528L. Screenshot 2024-12-01 at 11 16 04 am

I then tried to get clever and specify +over explicitly on the projection, but get the same result, missing everything west of the antimeridian: gdalwarp -of PNG -t_srs '+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 +over' -te -22094105.797 -11068715.659 -10018754.171 11068715.659 raw.tiff over.png Screenshot 2024-12-01 at 11 20 35 am

Trying with CENTER_LONG seems to have no effect gdalwarp -of PNG -t_srs 'EPSG:3857' --config CENTER_LONG 180 raw.tiff over.png Screenshot 2024-12-01 at 11 08 51 am

Unless I change the projection to EPSG:4326, in which case there's nothing east of the antimeridian?? gdalwarp -of PNG -t_srs 'EPSG:4326' --config CENTER_LONG 180 raw.tiff over.png Screenshot 2024-12-01 at 11 22 44 am

Versions and provenance

GDAL 3.10.0, released 2024/11/01 OSX via Brew

Additional context

No response

jratike80 commented 1 week ago

I think you must clip the raster into two halves.

gdalwarp -of PNG -t_srs EPSG:3857 raw.tiff out1.png -te_srs EPSG:4326 -180 -65 -90 70
gdalwarp -of PNG -t_srs EPSG:3857 raw.tiff out2.png -te_srs EPSG:4326  170 -65 180 70
Plantain commented 1 week ago

I think this should not have to be the case with CENTER_LONG or lon_wrap=180?

Because of some downstream decision making my tileserving speed is linearly related to the number of files I have to check.

jratike80 commented 1 week ago

I don't know how you could make such single wrapping image in the target CRS, EPSG:3857. The Proj library supports wrapping only with geographic coordinate systems, not with projected ones. See https://proj.org/en/stable/usage/projections.html

The +lon_wrap option can be used to provide an alternative means of doing longitude wrapping. It has only effect with operations that output angular coordinates, such as +proj=longlat

Plantain commented 1 week ago

Should either:


gdalwarp -of PNG -t_srs EPSG:3857 raw.tiff out.png -te_srs '+init=epsg:4326 +proj=longlat +ellps=WGS84
+datum=WGS84 +no_defs +towgs84=0,0,0 +over' -190 -65 -90 70

or


gdalwarp -of PNG -t_srs EPSG:3857 raw.tiff out.png -te_srs '+init=epsg:4326 +proj=longlat +ellps=WGS84
+datum=WGS84 +no_defs +towgs84=0,0,0 +lon_wrap=180' 170 -65 270 70

not work then?

Or is there any other way to achieve this? 391299766-fca1242e-0e07-4732-b747-17e10ab7894f

jratike80 commented 1 week ago

Please clarify what is your aim. Is it a single image like this, in EPSG:3857?

image
Plantain commented 1 week ago

Yes.

However, I also wanted to try centered on the antimeridian (like below) to compare. I was sure if they had would have the same computational complexity to generate or size to store.

Screenshot 2024-12-02 at 5 50 27 pm

jratike80 commented 1 week ago

The projected bounds of EPSG:3857 are

-20037508.34 -20048966.1
20037508.34 20048966.1

I do not believe that there is any supported way to centre EPSG:3857 at any other longitude than 0, or to go beyond the bounds of the valid coordinate range. So either you must split the image into two halves, or to select some other coordinate system and create a custom tiling schema for that CRS.

The two half-images can be merged into one with gdal_merge or gdalbuildvrt+gdal_translate if it is essential to have one image instead of two, but then you will have lots of nodata pixels in the middle.

Note: I have never done anything with a gestatationary CRS like this and I noticed that my own commands do not really work. The CRS is +proj=geos +sweep=x +lon_0=-137 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

jratike80 commented 6 days ago

Below are tested commands which creates two halves in EPSG:3857 which show correctly in QGIS. The images can then be converted into Web Mercator tiles according to the Google schema, if that was your final goal, The Google tiles are never split by the antimeridian so the end result will be the same from these two half-images of from a theoretical single image. Theoretical because I do not know if it is possible to create such.

gdalwarp -of GTiff -t_srs EPSG:3857 raw.tiff out1.tiff -co tiled=yes -co compress=LZW -te_srs EPSG:4326 -te -180 -65 -90 70
gdalwarp -of GTiff -t_srs EPSG:3857 raw.tiff out2.tiff -co tiled=yes -co compress=LZW -te_srs EPSG:4326 -te  170 -65 180 70
gdaladdo out1.tiff
gdaladdo out2.tiff

By creating Cloud Optimized GeoTIFFs the overviews and tiles are generated automatically but it is slower and I do not know if you would have any benefit from COG for your use case.

gdalwarp -of COG -t_srs EPSG:3857 raw.tiff out1_cog.tiff -co compress=LZW -te_srs EPSG:4326 -te -180 -65 -90 70
gdalwarp -of COG -t_srs EPSG:3857 raw.tiff out2_cog.tiff -co compress=LZW -te_srs EPSG:4326 -te  170 -65 180 70

I do not believe that you you have found any bug in GDAL so this issue may be closed soon. The right forum for asking questions about how to use GDAL is the gdal-dev mailing list. I know that in GitHub it is more convenient because it is easy to add screenshots etc.