Sunday, November 4, 2012

Zoomable image with Leaflet

Last week the European Southern Observatory (ESO) published a 9 giga-pixel image cataloging 84 million stars in the Central parts of the Milky Way. They even provided this image in a zoomable format at http://www.eso.org/public/images/eso1242a/zoomable/ .

It's really cool, but I'm not particularly fond of the Flash Plugin that they've used. It's a little bit slugish, and IMHO something like Google/Bing/Leaflet could provide a much better user experience.

That said, here's what I'll do: I'll get a smaller image from ESO, split it in tile images and display it on Leaflet.

To split the image I'm going to use an utility called GDAL2Tiles. I'll show two variants of this:
  • using the command-line tool
  • using a Graphical Interface called MapTiler
If you just want to keep things simple jump right ahead to the MapTiler version. I'm also showing the command-line version because it's more productive and suitable to batching scenarios.


I've obtained the image here. According to the description this is a "VST image of the star-forming region Messier 17".

1.a) Using GDAL2Tiles

Requirements:
  • Install Python2.7 (installer here)
  • Install GDAL-core (unofficial installer here)  (updated 12-02-2015, thanks Kalix)
  • Install GDAL-Python bindings (unofficial installer here)  (updated 12-02-2015, thanks Kalix)
Now, when running the GDAL command line, you have access to gdal2tiles.py.

You can check the help at the official page here.

Basically this is it.
gdal2tiles.py [-p profile] [-r resampling] [-s srs] [-z zoom]
              [-e] [-a nodata] [-v] [-h] [-k] [-n] [-u url]
              [-w webviewer] [-t title] [-c copyright]
              [-g googlekey] [-b bingkey] input_file [output_dir]
As you might suspect by the arguments, this utility is also able to automatically create a web viewer for the created tiles. In our case, as we want to do things manually and use Leaflet, we won't generate the viewers (hence, parameter 'w' will have the 'none' value).

Regarding the other parameters:
  • profile: the default is "mercator". In our case the image is not georeferenced, so the value should be "raster"
  • resampling: will use the default value "average"
  • SRS: Our image doesn't have any SRS, so no value here
  • resume mode ('e'): No, so don't use parameter 
  • Zoom level. See below:
The "zoom level" parameter is relatively easy to grasp:

On a Leaflet map (or Google/Bing), the furthest way zoom level (whole world on a single tile) is 0 (zero). This corresponds to a 256x256 square image. With each zoom level the square size is multiplied by 2, thus:
  • Zoom level 0: 256 px
  • Zoom level 1: 512 px
  • Zoom level 2: 1024 px
  • Zoom level 3: 2048 px
  • Zoom level 4: 4096 px
  • Zoom level 5: 8192 px
  • Zoom level 6: 16384 px
(and so on)

So, in my case the image is 16017x16017, so which zoom level to choose? Well, we can choose a max of 6, but this will make the map have 367 pixels of empty white space (16384-16017). So, we have three options here:
  1. Embrace the white bands and do nothing
  2. Resize the image so that it matches precisely a power of 2 (256, 512, ..., 16384)
  3. Convert it to a GeoTIFF file (which is georeferenced), where the corners match the world map corner coordinates.
For now I'll go for option 2, using another tool from GDAL to resize the image. So, to resize it to 16384 and keep it as a JPEG file I've used the following command:
gdal_translate -of JPEG -out size 16384 16384 eso1119a.jpg eso.jpg

meaning:
gdal_translate -of JPEG -out size «target witdth» «target height» «source file» «target file»

Now we're ready to run the GDAL2Tiles:
gdal2tiles.py -p raster -z 0-6 -w none eso.jpg

Just as simple as that. Wait a couple of minutes and you should have a folder with the following structure:
eso / «z» / «x» / «y».jpg



1.b) Using MapTiler

MapTiler has no requirements. Just download, run the installer and you're set to go.

I'll assume though that the input image is the one that was resized on the previous steps (eso.jpg).







 




You may notice that these settings map almost one-to-one to the GDAL2Tiles command-line version.

Anyway it has got everything inside the same package so you don't need to install Python and GDAL.

2) Using the tiles in Leaflet

Now for the really easy part, create the following html page:
<html>
<head>
    <link rel="stylesheet" 
          href="http://cdn.leafletjs.com/leaflet-0.6.4/leaflet.css" /> 
</head>
<body>
    <div id="map" style="width: 700px; height: 500px"></div>
 
    <script src="http://cdn.leafletjs.com/leaflet-0.6.4/leaflet.js"></script>

    <script>
        var map = L.map('map').setView([0, 0], 2);
        L.tileLayer('eso/{z}/{x}/{y}.jpg', {
            minZoom: 1,
            maxZoom: 6,
            attribution: 'ESO/INAF-VST/OmegaCAM',
            tms: true
        }).addTo(map);
    </script>
</body>
</html>

It's very important to set the "tms" parameter to "true". In previous Leaflet versions this was set with the "scheme" parameter.

3) End-Result

Check the web-page here.



34 comments:

  1. This comment has been removed by the author.

    ReplyDelete
  2. Thanks for this! I got this to work on OS X with the following modifications:
    -I used mac ports to install, with "sudo port install py-gdal"
    -I had to find where it was installed, it put it in:
    /opt/local/share/examples/py27-gdal/scripts/gdal2tiles.py
    -I had to run this for python to find py-gdal:
    sudo easy_install GDAL
    -It looks like -out size has been changed to -outsize, like this:
    gdal_translate -of JPEG -outsize 16384 16384

    ReplyDelete
  3. On ubuntu 12.04, i got the gda2tiles.py to work on the terminal. When i tried to run the command from PHP using exec() the following error was produced:

    python: /opt/lampp/lib/libz.so.1: no version information available (required by python)
    Traceback (most recent call last):
    File "/usr/bin/gdal2tiles.py", line 38, in
    from osgeo import gdal
    File "/usr/lib/python2.7/dist-packages/osgeo/__init__.py", line 21, in
    _gdal = swig_import_helper()
    File "/usr/lib/python2.7/dist-packages/osgeo/__init__.py", line 17, in swig_import_helper
    _mod = imp.load_module('_gdal', fp, pathname, description)
    ImportError: /usr/lib/i386-linux-gnu/libkrb5.so.26: undefined symbol: sqlite3_open_v2
    Error in sys.excepthook:
    Traceback (most recent call last):
    File "/usr/lib/python2.7/dist-packages/apport_python_hook.py", line 66, in apport_excepthook
    from apport.fileutils import likely_packaged, get_recent_crashes
    File "/usr/lib/python2.7/dist-packages/apport/__init__.py", line 1, in
    from apport.report import Report
    File "/usr/lib/python2.7/dist-packages/apport/report.py", line 20, in
    import apport.fileutils
    File "/usr/lib/python2.7/dist-packages/apport/fileutils.py", line 22, in
    from apport.packaging_impl import impl as packaging
    File "/usr/lib/python2.7/dist-packages/apport/packaging_impl.py", line 20, in
    import apt
    File "/usr/lib/python2.7/dist-packages/apt/__init__.py", line 21, in
    import apt_pkg
    ImportError: /usr/lib/i386-linux-gnu/libapt-pkg.so.4.12: undefined symbol: gzseek64

    Original exception was:
    Traceback (most recent call last):
    File "/usr/bin/gdal2tiles.py", line 38, in
    from osgeo import gdal
    File "/usr/lib/python2.7/dist-packages/osgeo/__init__.py", line 21, in
    _gdal = swig_import_helper()
    File "/usr/lib/python2.7/dist-packages/osgeo/__init__.py", line 17, in swig_import_helper
    _mod = imp.load_module('_gdal', fp, pathname, description)
    ImportError: /usr/lib/i386-linux-gnu/libkrb5.so.26: undefined symbol: sqlite3_open_v2

    ReplyDelete
  4. Very useful article for me.
    Thanks a lot!

    ReplyDelete
  5. Very informative and useful. Thank you.

    ReplyDelete
  6. If your leaflet page is greyed out check extension of the tiles. My tiles were .png's, but in the html provided it calls .jpg's, so change this line:
    L.tileLayer('eso/{z}/{x}/{y}.jpg', {
    to
    L.tileLayer('eso/{z}/{x}/{y}.png', {

    simple, stupid, but might be useful for somebody

    ReplyDelete
    Replies
    1. That was my problem lol
      Thank you so much

      Delete
  7. I am trying to build a leaflet based android map app that has to work offline. I have stored the tiles via MOBAC (Mobile Atlas Creator) onto my device's SD card (Internal Storage\tiles\mapnik\. I think my issue is the url for the tile layer. I have tried all combinations here ..L.tileLayer("/tiles/mapnik/"..Any suggestions on the url for referencing tiles on the device?

    Thanks!

    ReplyDelete
    Replies
    1. Hi Suma how did you download the map layer suing mobile atlas??

      Thnaks for your help

      Delete
    2. Using Mobile Atlas Creator I selected my region and zoom levels. It then creates a zip which which I extracted to the assets directory of my Android project. I then referenced the tiles as follows:

      L.tileLayer('file:///android_asset/map/Mapnik/{z}/{x}/{y}.png', {
      attribution: '© OpenStreetMap contributors'
      }).addTo(map);

      And the maps now show up offline!

      Delete
  8. Great article! Thank you. Minor correction: the `-out size` parameter in `gdal_translate` should be one word... e.g. `-outsize`

    ReplyDelete
  9. im stuck and need help, do you have time to help me ?

    ReplyDelete
  10. Thank you Pedro for a good tutoria!! I am developing a non-geographical map application on a page where all assets will come from a CDN which means no folderstructure and relativity in paths. Any idea how to hook Leaflet with that,Would it be possible to use an intermediate json object telling Leaflet the absolute urls of the tiles?

    ReplyDelete
  11. Thanks Pedro. I already have been tiling images using a map tiler that I built and I have been using the settings to work with Google Maps API. When I try to use your Leaflet viewer with my existing map tiles, my tile labeled 0_0_0 works fine (yes, I changed the slashes in the code to be underscores so it would match the Google naming structure), but when I zoom in to the next zoom level (which has 4 tiles), my bottom tiles are on top and my top tiles are on the bottom. Is there an easy way to configure this to work with google maps tiles naming structure?

    ReplyDelete
    Replies
    1. never mind, I figured out it was the 'tms' setting

      Delete
  12. Wow it's very instructive, thank you ! How can we integrate data to this map like the names of stars with zoom variation ?

    ReplyDelete
  13. Thats Great.. ! Thank you so much Pedro Sousa.. !
    But when i tried to create tiles using GDAL2Tiles using commands as you shown above, it is only result ouput like this " WindowsError: [Error 5] Access is denied: 'eso' "
    Is there any permission or something else to do this?
    Please help me.. Thanks..

    ReplyDelete
    Replies
    1. I'm sorry.. that was simple ..
      I just run it as administrator and boom...

      Delete
  14. For those wanting to recreate this exactly, the correct download links (as of Feb 2015) are now at:

    GDAL-core: download.gisinternals.com/sdk/downloads/release-1600-gdal-1-9-2-mapserver-6-2-0/gdal-19-1600-core.msi

    GDAL-Python: download.gisinternals.com/sdk/downloads/release-1600-gdal-1-9-2-mapserver-6-2-0/GDAL-1.9.2.win32-py2.7.msi

    ReplyDelete
  15. Thank you so much for this post!!!

    ReplyDelete
  16. this line might need to be improved:
    gdal_translate -of JPEG -out size 16384 16384 eso1119a.jpg eso.jpg
    to:
    gdal_translate -of JPEG -outsize 16384 16384 eso1119a.jpg eso.jpg

    Note the 'outsize' without a space. This is what worked for me. Thanks for this example!

    ReplyDelete
  17. Thankyou your tutorial are awesome. Is there any chance that I can put any marker to the map? Thank you for your answer.

    ReplyDelete
    Replies
    1. Did you find out how ?

      I am trying to get the coordinations with gimp and transform it to put the marker but it doesn't work ! :(

      Delete
  18. Awesome! At last! After trying out half a dozen different tutorials, here is (finally!) one that works — and apparently has been working for the past 4 years or so. Thank you so much for relieving my headaches :)

    @Ramadhan what I did in my case is simply to add:

    var popup = L.popup();
    function onMapClick(e) {
    popup
    .setLatLng(e.latlng)
    .setContent("You clicked the map at " + e.latlng.toString() + "\nZoom level is " + map.getZoom())
    .openOn(map);
    }
    map.on('click', onMapClick);


    That way, I can just point and click to where I wish to place a marker, copy & paste the lat/lng, and that's it :)

    ReplyDelete
  19. Thanks for this! Everything worked perfectly.

    ReplyDelete
  20. Very nice, thanks!
    Just to mention that I've found libvips very useful for tile generation (https://github.com/libvips/libvips)

    ReplyDelete