El problema es que el array `lunar` no es un array de dos dimenciones porque la imagen no es una imagen en monocromática, sino que es en color. Esto hace que cada pixel deba representarse mediante un array con el valor de cada canal, por lo que tiene  tres dimensiones, por ejemplo una imágen de 2 x 2 RGBA/BGRA sería algo como:

    [[[0.2, 0.3, 1.0, 1.0],
      [0.9, 1.0, 0.6, 1.0]],
     [[1.0, 1.0, 0.4, 1.0],
      [1.0, 1.0, 1.0, 0.9]]]

en la que cada pixel es un array de la forma `[R, G, B, A]`/`[B, G, R, A]`.

Por lo tanto  `shape` retorna una tupla de tres elementos, el primero son las filas, el segundo las columnas y el tercero son el valor para cada canal del pixel:

Debes hacer por tanto:

    M, N, CHANNELS = lunar.shape

o si no quieres para nada el número de canales:

    M, N, _ = lunar.shape

o una forma general para solo desempaquetar los dos primeros items y descartar el resto sin importar el número: 

    M, N, *_ = lunar.shape

Obviamente esto no solventa el problema si necesitas una imagen monocromática, en tal caso puedes convertirla combinarla manualmente con NumPy ([ITU-R 601-2][1]):

    lunar_gray = np.dot(lunar[...,: 3], [299/1000, 587/1000, 114/1000])

 
O usa directamente `pillow`.


  [1]: https://en.wikipedia.org/wiki/Rec._601