1

I working on a project need to tackle with deal with geotiff image file. I'm new to leaflet and geotiff image.I got an image from client and i need to overlay it to leaflet map from geotiff. I used

gdal2tiles.py -p raster [image.tif]
to generater tiles map. And i tried to display it on leafletmap, however I got this displayDisplaying from react-leaflet

I wondering anyone have idea to figure out, how to find out the optimal zoom level for convertion from geotiff to raster files without lossing it's quality?Here is my gdalinfo from my geotiff file.

Driver: GTiff/GeoTIFF
Files: Geotiff.tif
Size is 25070, 18673
Coordinate System is:
PROJCS["WGS 84 / UTM zone 54N",
    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"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",141],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32654"]]
Origin = (297884.090830000001006,3955249.192560000345111)
Pixel Size = (0.071710000000000,-0.071710000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=pix4dmapper
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  297884.091, 3955249.193) (138d45'55.76"E, 35d43'13.77"N)
Lower Left  (  297884.091, 3953910.152) (138d45'56.97"E, 35d42'30.34"N)
Upper Right (  299681.861, 3955249.193) (138d47' 7.26"E, 35d43'15.09"N)
Lower Right (  299681.861, 3953910.152) (138d47' 8.46"E, 35d42'31.66"N)
Center      (  298782.976, 3954579.672) (138d46'32.11"E, 35d42'52.71"N)
Band 1 Block=25070x1 Type=Float32, ColorInterp=Gray
  NoData Value=-10000
  • I can't say about zoom levels but it seems that your Leaflet map is configured to have tile origin in wrong place. In TMS the origin is at the lower left corner wiki.osgeo.org/wiki/Tile_Map_Service_Specification. – user30184 Jan 16 '18 at 15:07
  • Thank you for reply @user30184. I think you are correct. The gdal2tiles.py install by default , it can't generate a TMS for leaflet. I found a modified version of gdal2tiles in github it perfectly shows how to find the optimal zoom level for converting to TMS. I managed to display my geotiff in leaflet, but i still got multiple display when I zoomed to high level and a corrupted TMS. Do you know what is the problem with this? I just wondering is the corrupted image caused by the compression LZW – Chee-siang Leow Jan 18 '18 at 8:12
  • Is the Leaflet map that you talk about the one that gdal2tiles creates with the -w leafletoption? – user30184 Jan 18 '18 at 9:50
  • yes, however , my gdal2tiles.py install by default the option Web viewer options: Options for generated HTML viewers a la Google Maps -w WEBVIEWER, --webviewer=WEBVIEWER Web viewer to generate (all,google,openlayers,none) - is show as this. I can't choose leaflet . Is this relate with the version of gdal? – Chee-siang Leow Jan 18 '18 at 15:36
  • I suppose so. But how can you use the Leaflet option if you do not have it available? – user30184 Jan 18 '18 at 17:24
1

I had a case like this, went through gdal2tiles.py source for inspiration and wrote this Python script:

from osgeo import gdal

def get_optimal_zoom_level (geo_tiff_file_name):
  geo_tiff = gdal.Open(geo_tiff_file_name)
  geo_transform = geo_tiff.GetGeoTransform()
  degrees_per_pixel = geo_transform[1]
  radius = 6378137
  equator = 2 * math.pi * radius
  meters_per_degree = equator / 360
  resolution = degrees_per_pixel * meters_per_degree
  pixels_per_tile = 256
  zoom_level = math.log((equator/pixels_per_tile)/resolution, 2)
  MAX_ZOOM_LEVEL = 20
  optimal_zoom_level = min(math.floor(zoom_level), MAX_ZOOM_LEVEL)
  return optimal_zoom_level

In your case you'd call

get_optimal_zoom_level('[image.tif]') 

and it should return deepest zoom level that still makes sense to tile at.

Note: This is meant to work for imagery not too far from the Equator and WGS 84 coordinate system. You may want to tweak it if your imagery is closer to poles.

Edit: My GeoTIFF unit was degree, so this code assumes unit is a degree. Your unit is meter, you'll have to tweak it for that, too.

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.