mapnik / mapnik

Mapnik is an open source toolkit for developing mapping applications
http://mapnik.org
GNU Lesser General Public License v2.1
3.64k stars 827 forks source link

Resampling broken with gdal raster reprojection #1034

Closed tomtaylor closed 10 years ago

tomtaylor commented 12 years ago

Resampling doesn't work correctly for gdal images which are being reprojected on the fly. The work around for now is to use gdalwarp to adjust them to your map's projection before rendering.

woodpeck commented 12 years ago

I was about to open a new issue when I found this old one which has perhaps been sidelined because it is not very specific. I am fighting raster reprojection on a regular basis and here's a detailed example, based on SRTM hillshading.

I have the following style that simply displays hillshading and nothing else:

<Map background-color="#fff" 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 +over +no_defs">
<Style name="hillshade">
   <Rule>
      <RasterSymbolizer opacity="0.4" scaling="bilinear8" mode="normal" />
   </Rule>  
</Style>

<!-- SRS identical to above -->
<Layer name="srtm" status="on" 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 +over +no_defs">
   <StyleName>hillshade</StyleName>
   <Datasource>
      <Parameter name="type">gdal</Parameter>
      <Parameter name="file">hillshade.tiff</Parameter>
      <Parameter name="format">tiff</Parameter>
   </Datasource>
</Layer>
</Map>

The source image for this test is a 32x32 image cropped from SRTM hillshading, shown magnified to 256x256 here:

raw hillshading

The desired result, and indeed the result with the current latest 2.1.0+dev20120816 build, is this:

correct result

However, when I add a blank to the SRS definition in the <Layer> tag and thereby force Mapnik to do a raster reprojection (even if the mathematical outcome of that reprojection is a null operation), I get this result with ugly white speckles:

bad result

(NB, this is an improvement over 2.1.0+dev20120703 which skipped smoothing altogether in this situation and created a very blocky output similar to the raw TIFF posted at the top.)

My issue with this is not that Mapnik runs the projection code at all, not seeing that both projections are indeed the same; I chose two different-looking-but-identical projections for testing only. My issue is that I actually want to use Mapnik in a situation where raster reprojection is not avoidable (namely, a WMS server that suports dozens of projections and it would be a major hassle to pre-project every raster into every possible projection).

I have prepared a little manual test case consisting of a shell script that you can run to produce the above "correct" and "broken" results with nik2img, here: http://www.remote.org/frederik/tmp/hillshadeprojection.zip - would be interesting to hear if others can confirm the problem or maybe even have a workaround (that is different from "just use gdalwarp...").

springmeyer commented 12 years ago

Hi @woodpeck - yes, we are aware that raster reproject is extremely buggy atm. Could you please create a new issue so this does not get lost. Or just fix the image links in this one.

woodpeck commented 12 years ago

Oops, sorry, fixed the image links!

woodpeck commented 11 years ago

Just tried this again today to see if anything had accidentally changed in the mean time. While the "broken" image is unchanged, the one that should be correct now looks like this (in 2.2.0+dev20130605.git.049b7b2):

broken even without reprojecting

springmeyer commented 11 years ago

hmm, odd spots, thanks for reporting.

woodpeck commented 11 years ago

I have investigated this further and found that the spots you're seeing are "nodata" pixels, and the problem goes away if you fix the nodata values in the source. So there is a bug but it seems to apply to the nodata masking and not reprojection in general. - Out of interest, I see that Mapnik has its own reprojection code for raster sources which gets applied even when the source is GDAL which also offers (potentially more sophisticated?) reprojection. Has using GDAL reprojection been attempted and rejected because of speed/quality, or has it never been tried?

springmeyer commented 11 years ago

Oh these are nodata pixels: interesting since that may point to a very specific bug. The warping code both warps and scales. I recently fixed a few alpha handling bugs in the scaling code but those may not have been ported to the separate scaling code in the warper. Is your testcase link above still working? Iif my guess is right then it should be a quick fix when I find time to look on Monday.

As far as why we warp after gdal: I think @albertov was working on adding overview reading support to the gdal plugin at the same time he was working on reprojection and ran into issues getting both working together. Design wise it also made sense to warp outside the plugin because plugins currently do not actually know the target map srs to be able to set up a transform. That said given that our custom warping code still is not production ready I do think it would be worthwhile to revisit trying to leverage gdal for warping.

woodpeck commented 11 years ago

Yes, my test case is still "working" in that the hillshade.tiff supplied with the test case does contain nodata pixels which currently produce spots. The test case used to produce different results depending on whether source and target projections were "only mathematically identical" or "even literally identical", but that distinction seems to have vanished at some point; run against trunk both results show the same spots.

When the source hillshade.tiff is viewed with software that is not "nodata aware" and ignores the TIFF tag, it looks normal (as the "nodata value" is 181 which is a light grey). When viewed with a "nodata aware" viewer, the nodata pixels show up in some implementation-dependent way.

Andrey-VI commented 11 years ago

Just faced with the same issue. I have no reprojection, but I have scaling option enabled on some levels. And on that levels I see same bug related to "NODATA" value (as seen on the left side of the image below): _scaling_bug

For the left side I used the clipped hillshade layer and non-clipped for the right.

springmeyer commented 11 years ago

@woodpeck - so I downloaded and can replicate. However, my sense is that an image with white spots is the correct behavior. This is because those are, like you notice, no data pixels, and no data pixels are considered transparent in Mapnik: given your background is white then the white should shine through. The reason you are getting gray spots is because you are using bilinear8 which is triggers custom scaling code and only exists for back compatibility. bilinear8 is not handling pre-multipled alpha correctly. We could consider fixing it, but ideally you could just move to using bilinear for scaling instead. Is there a specific reason you are using bilinear8?

dieterdreist commented 11 years ago

I also encountered something that might be related. On OSX rendering tiles with mapnik 2.2.0 (built yesterday with brew install mapnik) I get segfaults when rendering tiles from gdal-tifs with generate_tiles.py every 5-30 tiles (it happens on an irregular basis) while using bilinear8, but it doesn't happen when I switch to bilinear. Unfortunately visually I prefer the bilinear8 output (my tifs contain dense hatches and with bilinear scaling I get more undesired artefacts). Btw.: the same issue was already in mapnik2.1 (segfaults with bilinear8).

Andrey-VI commented 10 years ago

@albertov , @springmeyer , any progress on this?

Andrey-VI commented 10 years ago

Small testcase here: http://91.208.39.38/~andrey/tiles/temp/tc_nodata.zip With scaling="near": nodata-near

With any other scaling option (scaling="quadric" here): nodata-quadric

Similar to https://github.com/mapnik/mapnik/issues/1508

springmeyer commented 10 years ago

@Andrey-VI - this issue does look similiar to #1508, but is I'm fairly sure is rather due to the comp-op:grain-merge you are using. So I have created https://github.com/mapnik/mapnik/issues/2067 to track this.

springmeyer commented 10 years ago

@dieterdreist - in reply to your comment above, I think you are hitting #1529. Can you please provide images there that help us understand exactly how using bilinear8 results in better looking maps? I'd prefer just to remove bilinear8.

springmeyer commented 10 years ago

@woodpeck - I'm going to close this issue. If you feel like my analysis in this comment is incorrect please re-open a new issue that clarifies what we should fix. Thanks!

btw, raster reprojection bugs still exist, but let's discuss them on other tickets like #1049

Andrey-VI commented 10 years ago

this issue does look similiar to #1508, but is I'm fairly sure is rather due to the comp-op:grain-merge you are using.

No, it's not due to comp-op. I've changed the testcase style to this:

<Map background-color="#faf9f5" srs="+proj=merc +datum=WGS84 +over" maximum-extent="9717000, 6254000, 9742000, 6262000" minimum-version="2.2.0">
  <Style name="hillshade">
    <Rule>
      <RasterSymbolizer scaling="quadric"/>
    </Rule>
  </Style>

  <Layer name="hillshade_zl11" status="on" srs="+proj=merc +datum=WGS84 +over">
    <StyleName>hillshade</StyleName>
    <Datasource>
      <Parameter name='type'>gdal</Parameter>
      <Parameter name='file'>hs_zl11_23_48.tif</Parameter>
    </Datasource>
  </Layer>
</Map>

and the bug is still here…

dieterdreist commented 10 years ago

sorry it took me some time to reply. I finally discovered that there were lots of nice resampling algorithms and tested all of them. The best results across all zoomlevels I got with sinc, so it didn't matter any more what was the bilinear problem. I am now recreating tiles with bilinear and bilinear8 to show you the difference. (Update: sorry, had a problem with my mapnik supposedly because of the upgrade to Maverick, am still struggling to get it back working)

springmeyer commented 10 years ago

@dieterdreist - okay, in the meantime I plan to remove bilinear8: #2076

dieterdreist commented 10 years ago

bilinear-bilinear8

@springmeyer I was finally able to recreate some tiles in bilinear and bilinear8, see attachment. Generally the bilinear stuff looks better, but in some situations and with a lot of scale reduction the bilinear8 seem to be sharper and of better contrast, while the bilinear seems a bit "blurred". Personally I am happy since I discovered the various other resampling options (sinc etc.)

springmeyer commented 10 years ago

@dieterdreist - interesting finding regarding sharpness. You may be able to achieve this with bilinear by applying a sharpen image-filter?

<Style image-filters="sharpen()" >
...
</Style>
``