3

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:

Original

Projected:

Projected

Both:

Both together

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

enter image description here

-tr 30 30

enter image description here

-tr 8 8

enter image description here

-tr 7.0832868662 7.0832868662

enter image description here

  • the output file is bigger because you haven't used any compression -co COMPRESS=LZW should help – Ian Turton Feb 15 '17 at 14:09
  • I can use LZW however it reduces the performance of geoserver's tile rendering, that's why I haven't used that. The reason I am reprojecting the data to another projection is to improve Geoserver's rendering time. – nilgun Feb 15 '17 at 14:12
  • 1
    for a 30gb file compression will help more than it hinders, also add overviews – Ian Turton Feb 15 '17 at 14:15
  • 1
    Try also gdalwarp -t_srs EPSG:3857 -co 'TILED=YES' -co 'BLOCKXSIZE=256' -co 'BLOCKYSIZE=256' original.tif projected.tifwithout target size. See if pixels are shifted and tell us what is the pixel size when GDAL can decide it. – user30184 Feb 15 '17 at 14:22
  • 1
    @nilgun I don't know what QGIS does, but GeoServer does it fast because most of the time the transformation can be assimilated to a simple scale (with different factors in x and y) without visible loss of accuracy. The complete trigonometric math is used only if the area being displayed is very large (low zoom levels) – Andrea Aime Feb 18 '17 at 15:27
2
+50

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.

  • Original pixel size is (0.000395841574135,-0.000395841574135). I am using a linestring from the center of the tif to calculate size of original pixels in 3857 as follows: SELECT ST_Length( ST_Transform( ST_GeomFromEWKT('SRID=4326;LINESTRING(35.2433229 38.9565576,35.2433229 38.9569534416)'), 3857 ) );. This yields 56.6662949295714. Since you suggested multiples, I divided this number by 8 and used -tr 7.0832868662 7.0832868662. The result is not as close (to the original 4326 on the fly projected result) as -tr 8 8. – nilgun Feb 16 '17 at 14:36
  • Do I need to manually check all my other files (there are like 30 more) to find out 8 8 works ok? – nilgun Feb 16 '17 at 14:36
  • I am already using GeoWebCache. But for high zoom levels (where it does not hit the cache) it takes like 18-25 seconds for geoserver to create tiles from my rasters with Type:Float32 with native:4326 srs, requested with srs:3857 from leaflet. – nilgun Feb 16 '17 at 14:40
  • If the size is around 56, I was thinking you could use, say, 20 (around 1/3) or maybe 30. The problem, I think, is that with a cell size of 56, and cells not lining up perfectly, your max error is pretty big. You can lower it by using smaller cells, but perhaps 8 is too small... – Alex Leith Feb 16 '17 at 22:49
  • 1
    Hey @nilgun, I reckon that GeoServer is rendering these at a very high resolution, compared with what you're doing. The pixels that GeoServer renders at end up being the same as your -tr flag, so if it's doing a very zoomed in image, the -tr might be 0.01 m, say, and in that case, it'll line up very well! – Alex Leith Feb 21 '17 at 23:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.