gdalwarp tiff reprojection from EPSG:4326 to EPSG:3857

Question

I have a geoTiff file with the following gdalinfo output:

Driver: GTiff/GeoTIFF Files: original.tif Size is 48377, 15906 Coordinate System is: GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY["EPSG","4326"]] Origin = (25.668509000000000,42.104685623790893) Pixel Size = (0.000395841574135,-0.000395841574135) Metadata: AREA_OR_POINT=Area DataType=Generic Image Structure Metadata: COMPRESSION=LZW INTERLEAVE=BAND Corner Coordinates: Upper Left ( 25.6685090, 42.1046856) ( 25d40' 6.63"E, 42d 6'16.87"N) Lower Left ( 25.6685090, 35.8084295) ( 25d40' 6.63"E, 35d48'30.35"N) Upper Right ( 44.8181368, 42.1046856) ( 44d49' 5.29"E, 42d 6'16.87"N) Lower Right ( 44.8181368, 35.8084295) ( 44d49' 5.29"E, 35d48'30.35"N) Center ( 35.2433229, 38.9565576) ( 35d14'35.96"E, 38d57'23.61"N) Band 1 Block=128x128 Type=Byte, ColorInterp=Gray NoData Value=15 Image Structure Metadata: NBITS=4

I want to reproject it to EPSG:3857 for native (fast) rendering on Geoserver. This is the command I use to reproject it:

gdalwarp -t_srs EPSG:3857 -co 'TILED=YES' -co 'BLOCKXSIZE=256' -co 'BLOCKYSIZE=256' -ts 48377 15906 original.tif projected.tif

Below, you can see the original, projected and both tiffs on top of each other:

Original:

Projected:

Both:

As you can see the projected raster is not aligned with the original one.

I have read that tapcan be used to align pixels, and tried this:

gdalwarp -t_srs EPSG:3857 -co 'TILED=YES' -co 'BLOCKXSIZE=256' -co 'BLOCKYSIZE=256' -tap -tr 8 8 original.tif projected.tif

Now this outputs a tiff file with a very little shift (about 2m) which I am ok with. But the processing is very slow and the resulting file is 30GB.

Giving large tr values such as -tr 40 40, fixes the slow processing problem, but the output tiles are no longer have the same value as the original ones because of the resampling I guess.

I also want to note that, tiles from the above picture are requested using EPSG:3857 srs from Geoserver. Geoserver can project native EPSG:4326 tiles to EPSG:3857 on the fly with accurate precision, but I can not produce the same result with gdal.

EDIT

These are some screenshots for different target resolutions:

-tr 20 20

-tr 30 30

-tr 8 8

-tr 7.0832868662 7.0832868662

Answer 1

What you could do is measure the size of your original pixels in your target SRS, and make your target resolution a multiple of those. (You'll need to pick something that works for most of the area, as it will change at high latitudes).

Otherwise, your testing with 8 x 8 sounds like it worked ok, so experiment with compression to get the file size right.

And finally, you could just rely on GeoServer and use GeoWebCache to convert it for you. This will result in something fast, as is what you want to be using anyway in production. So, even if you do optimise your source data, allowing GWC to pre-cache or cache on demand should make it as good as possible.

EDIT:

What's going on here is that your TARGET raster is not going to be aligned with your SOURCE raster. As you note, using a low resolution target means that the error in alignment is small, but it't not perfect. This is a fundamental fact of the two grids you're making, in two quite different coordinate reference systems, and unless you make your grid very small, you'll always have errors. (And unless that small grid lines up perfectly with every cell of the source, it'll still have errors, really, they'll just be small...)

Something to ask is: how accurate is your data. You're mapping at a ~50 m grid, recall, so it's not going to be very accurate, and is certainly not all that precise. A target grid that lines up within maybe 10%, or 5 m, could be acceptable. But that's up to you.

Looking at your images, the 20 m and 30 m grids line up good enough, in my view. It's not a systematic error, i.e., you've got errors in all different directions.