How to use Leaflet to display a map in Lambert Conformal Conic without warping to EPSG:3857

Question

I would like to display a map in leaflet which is in Lambert Conformal Conic. The CRS of this map is: '+proj=lcc +lat_0=38 +lon_0=-100 +lat_1=45 +lat_2=33 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'

I do not want to warp the image into EPSG:3857 as I want the image on the web to be identical to the image as it would be printed in real life.

I've been playing around with using proj4leaflet, but I am struggling trying to figure out the resolution and origin. I don't know if this is my problem, but my pixel size is different for the x and y-direction.

Pixel Size = (262.477506424478293,-262.465501517241137)

My original image size is:

Size is 18509, 11441

So, given a 262.477506 and -262.4655 I would think at zoom level 0 with a 256x256 tile, my resolution (projection units (meters) to pixels) would be the total width in meters 18509*262.477506/256 would be 18,977.3287443516 meters to a pixel. But that just isn't working?

My steps are to:

gdal2tiles.py --profile raster U_S_VFR_Wall_Planning_Chart.tif -z 0-5 --xyz

Then I am using this:

var crs = new L.Proj.CRS("",'+proj=lcc +lat_0=38 +lon_0=-100 +lat_1=45 +lat_2=33 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs',
  {
    origin: ??
    resolutions: ??
  }) 

    var mymap = L.map('mapid', { 
    minZoom: 0, 
    maxZoom: 5,
    continuousWorld: true,
    crs: crs, 
    tms: true,
    noWrap: true
}).setView([26,-80], 0);
    L.tileLayer('./U_S_VFR_Wall_Planning_Chart/{z}/{x}/{y}.png', {
    }).addTo(mymap);

When looking at the tiles they appear to be correct:

Answer 1

I got this all worked out. However, I still do not quite understand how exactly it worked with respect to the resolutions.

The origin can be taken right from the gdalinfo output.

Origin = (-2189360.538664901163429,1531470.088855510344729)

It's the lower left corner to the upper right corner.

For the resolutions, I looked at how many zoom levels were created and then worked backwards having the final zoom layer been the actual pixel size from the gdalinfo. Each higher zoom level is just double that until you get to zero.

ex: For our pixel size of ~262.477, each zoom level will double that.

[33597.1208223332214,16798.5604111666107,8399.28020558330536,4199.64010279165268,2099.82005139582634,1049.91002569791317,524.955012848956585,262.477506424478293],

What I don't quite fully understand is how gdal2tiles figures out what the max zoom level is. I see that the highest zoom level has the resolution of the pixel width, why is it zoom level 7 and not 8?

Basically all this information is in the openlayers.html file generated by gdal2tiles. I just copied it into my script.

Here is the final solution:

  var crs = new L.Proj.CRS("", '+proj=lcc +lat_0=38 +lon_0=-100 +lat_1=45 +lat_2=33 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs',
      {
        origin: [-2071376.899753200123087,1348662.866750749992207],
        resolutions: [ 8399.28020479999941,
4199.64010239999971,
2099.82005119999985,
1049.91002559999993,
524.955012799999963,
262.477506399999982 
        ]
      })

    var mymap = L.map('map', {
      minZoom: 0,
      maxZoom: 5,
      crs: crs,
    }).setView([26, -80], 0);
    L.tileLayer('./tiles/{z}/{x}/{y}.png', {
    }).addTo(mymap);