How to make smoother 2d maps out of a huge (Mbs-scale) .txt file?

Question

I have a .txt file that contains five columns of data in which columns 1 and 2 correspond to projected distance along dimension 1 and 2 (pixel scale of 0.50) and then columns 3, 4 and 5 are values of an arbitrary (otherwise unknown) function in projections 0 "z", 1 "x" and 2 "y" respectively. (Original file is ~1e6 rows long but here is a small version of it that I managed to upload: m11q.txt) I would like to transform this .txt table into actual 2d maps in which the color strength of each pixel would correspond to the magnitude of the function in that pixel. This means that I would have three planes xy, yz and xz in which each pixel is filled with a color.

My values for columns 3,4 and 5 are mostly zero with few of them on the order of $10^{14}$ to $10^{21}$. I am interested in having three contour levels. Let's say if the function value is smaller than $10^{14}$, the color should be black, if the function value is between $10^{14}$ and $10^{15.5}$, the color should be yellow, and if the function value is between $10^{15.5}$ and $10^{17}$, the color should be orange, if the function value is between $10^{17}$ and $10^{20.3}$, the color should be red and finally if above $10^{20.3}$ the color should be blue. In other words, my 2d map is in fact a set of 5 contour plots that are produced out of .txt table. Here is my original version of the code:

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate
import matplotlib.figure as fig
import matplotlib.cm as cm
import matplotlib.colors as colors
import matplotlib.colorbar as cb
import matplotlib.ticker as tk 
import pylab as pl

x, y, covering_fraction_z, covering_fraction_x, covering_fraction_y = np.genfromtxt(r'm11q.txt', unpack=True)


nInterp = 4000
xi, yi = np.linspace(x.min(), x.max(), nInterp), np.linspace(y.min(), y.max(), nInterp)
xi, yi = np.meshgrid(xi, yi)

z0 = scipy.interpolate.griddata((x, y), covering_fraction_z, (xi, yi), method='nearest')
z1 = scipy.interpolate.griddata((x, y), covering_fraction_x, (xi, yi), method='nearest')
z2 = scipy.interpolate.griddata((x, y), covering_fraction_y, (xi, yi), method='nearest')

LAF_covering_fraction_XY = np.count_nonzero((np.power(10, 14) <= z0) & (z0 < np.power(10, 17))) * 1.0 / z0.size
LAF_covering_fraction_YZ = np.count_nonzero((np.power(10, 14) <= z1) & (z1 < np.power(10, 17))) * 1.0 / z1.size
LAF_covering_fraction_XZ = np.count_nonzero((np.power(10, 14) <= z2) & (z2 < np.power(10, 17))) * 1.0 / z2.size
LAF = np.sqrt(np.power(LAF_covering_fraction_XY, 2)+np.power(LAF_covering_fraction_YZ, 2)+np.power(LAF_covering_fraction_XZ, 2))/np.sqrt(3)
LAF_sigma = max([abs(LAF-LAF_covering_fraction_XY),abs(LAF-LAF_covering_fraction_YZ),abs(LAF-LAF_covering_fraction_XZ)])
LLS_covering_fraction_XY = np.count_nonzero((np.power(10, 17) <= z0) & (z0 < np.power(10, 20.3))) * 1.0 / z0.size
LLS_covering_fraction_YZ = np.count_nonzero((np.power(10, 17) <= z1) & (z1 < np.power(10, 20.3))) * 1.0 / z1.size
LLS_covering_fraction_XZ = np.count_nonzero((np.power(10, 17) <= z2) & (z2 < np.power(10, 20.3))) * 1.0 / z2.size
LLS = np.sqrt(np.power(LLS_covering_fraction_XY, 2)+np.power(LLS_covering_fraction_YZ, 2)+np.power(LLS_covering_fraction_XZ, 2))/np.sqrt(3)
LLS_sigma = max([abs(LLS-LLS_covering_fraction_XY),abs(LLS-LLS_covering_fraction_YZ),abs(LLS-LLS_covering_fraction_XZ)])
DLA_covering_fraction_XY = np.count_nonzero(np.power(10, 20.3) <= z0) * 1.0 / z0.size
DLA_covering_fraction_YZ = np.count_nonzero(np.power(10, 20.3) <= z1) * 1.0 / z1.size
DLA_covering_fraction_XZ = np.count_nonzero(np.power(10, 20.3) <= z2) * 1.0 / z2.size
DLA = np.sqrt(np.power(DLA_covering_fraction_XY, 2)+np.power(DLA_covering_fraction_YZ, 2)+np.power(DLA_covering_fraction_XZ, 2))/np.sqrt(3)
DLA_sigma = max([abs(DLA-DLA_covering_fraction_XY),abs(DLA-DLA_covering_fraction_YZ),abs(DLA-DLA_covering_fraction_XZ)])


z0plot = plt.subplot(221)
plt.imshow(z0, vmin=0.0001+covering_fraction_z.min(), vmax=covering_fraction_z.max(), origin='lower',
           extent=[x.min(), x.max(), y.min(), y.max()],
                norm=colors.LogNorm())
plt.ylabel("Y [kpc]")
z0plot.set(aspect=1)


axparameters = plt.subplot(222)
frame = plt.gca()
frame.axes.get_xaxis().set_ticks([])
frame.axes.get_yaxis().set_ticks([])
axparameters.set(aspect=1)
plt.text(0.02, 0.2, 'LAF=%.5f$\pm$%.5f \n\n\n LLS=%.5f$\pm$%.5f \n\n\n DLA=%.5f$\pm$%.5f' % (LAF, LAF_sigma, LLS, LLS_sigma, DLA, DLA_sigma))
plt.colorbar()
plt.tight_layout()

z1plot = plt.subplot(223)
plt.imshow(z1, vmin=0.0001+covering_fraction_y.min(), vmax=covering_fraction_y.max(), origin='lower',
           extent=[x.min(), x.max(), y.min(), y.max()],
                norm=colors.LogNorm())
plt.xlabel("X [kpc]")
plt.ylabel("Z [kpc]") 
z1plot.set(aspect=1)



z2plot = plt.subplot(224)
plt.imshow(z2, vmin=0.0001+covering_fraction_x.min(), vmax=covering_fraction_x.max(), origin='lower',
           extent=[x.min(), x.max(), y.min(), y.max()],
                norm=colors.LogNorm())
plt.xlabel("Y [kpc]")
z2plot.set(aspect=1)


myfig = fig.Figure(figsize=(8, 8), dpi=900, facecolor='w', edgecolor='k', tight_layout=True)
plt.suptitle('HI Column Density Maps', size=20)
plt.tight_layout(rect=[0, 0, 1, 0.9])
plt.savefig('m11q.png', dpi=900)
plt.show()

Which produces the following map

---

Here is the code written by Xevaquor as an answer to my original post:

# -*- coding: utf-8 -*-    
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate
#import matplotlib.figure as fig
import matplotlib.cm as cm
import matplotlib.colors as colors
import matplotlib.colorbar as cb
import matplotlib.ticker as tk 
import pylab as pl
import matplotlib.ticker as ticker
from matplotlib import gridspec

x, y, covering_fraction_z, covering_fraction_x, covering_fraction_y = np.genfromtxt(r'm11q.txt', unpack=True)

def fmt(x, pos):
    return '10e{}'.format(np.log10(x))

nInterp = 4300
xi, yi = np.linspace(x.min(), x.max(), nInterp), np.linspace(y.min(), y.max(), nInterp)
xi, yi = np.meshgrid(xi, yi)

z0 = scipy.interpolate.griddata((x, y), covering_fraction_z, (xi, yi), method='nearest')
z1 = scipy.interpolate.griddata((x, y), covering_fraction_x, (xi, yi), method='nearest')
z2 = scipy.interpolate.griddata((x, y), covering_fraction_y, (xi, yi), method='nearest')

LAF_covering_fraction_XY = np.count_nonzero((np.power(10, 14) <= z0) & (z0 < np.power(10, 17))) * 1.0 / z0.size
LAF_covering_fraction_YZ = np.count_nonzero((np.power(10, 14) <= z1) & (z1 < np.power(10, 17))) * 1.0 / z1.size
LAF_covering_fraction_XZ = np.count_nonzero((np.power(10, 14) <= z2) & (z2 < np.power(10, 17))) * 1.0 / z2.size
LAF = np.sqrt(np.power(LAF_covering_fraction_XY, 2)+np.power(LAF_covering_fraction_YZ, 2)+np.power(LAF_covering_fraction_XZ, 2))/np.sqrt(3)
LAF_sigma = max([abs(LAF-LAF_covering_fraction_XY),abs(LAF-LAF_covering_fraction_YZ),abs(LAF-LAF_covering_fraction_XZ)])
LLS_covering_fraction_XY = np.count_nonzero((np.power(10, 17) <= z0) & (z0 < np.power(10, 20.3))) * 1.0 / z0.size
LLS_covering_fraction_YZ = np.count_nonzero((np.power(10, 17) <= z1) & (z1 < np.power(10, 20.3))) * 1.0 / z1.size
LLS_covering_fraction_XZ = np.count_nonzero((np.power(10, 17) <= z2) & (z2 < np.power(10, 20.3))) * 1.0 / z2.size
LLS = np.sqrt(np.power(LLS_covering_fraction_XY, 2)+np.power(LLS_covering_fraction_YZ, 2)+np.power(LLS_covering_fraction_XZ, 2))/np.sqrt(3)
LLS_sigma = max([abs(LLS-LLS_covering_fraction_XY),abs(LLS-LLS_covering_fraction_YZ),abs(LLS-LLS_covering_fraction_XZ)])
DLA_covering_fraction_XY = np.count_nonzero(np.power(10, 20.3) <= z0) * 1.0 / z0.size
DLA_covering_fraction_YZ = np.count_nonzero(np.power(10, 20.3) <= z1) * 1.0 / z1.size
DLA_covering_fraction_XZ = np.count_nonzero(np.power(10, 20.3) <= z2) * 1.0 / z2.size
DLA = np.sqrt(np.power(DLA_covering_fraction_XY, 2)+np.power(DLA_covering_fraction_YZ, 2)+np.power(DLA_covering_fraction_XZ, 2))/np.sqrt(3)
DLA_sigma = max([abs(DLA-DLA_covering_fraction_XY),abs(DLA-DLA_covering_fraction_YZ),abs(DLA-DLA_covering_fraction_XZ)])


fig = plt.figure(figsize=(8, 8))

gs = gridspec.GridSpec(2, 3, width_ratios=[5, 5, 1], height_ratios=[1,1,1]) 
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[3])
ax2 = plt.subplot(gs[4])

#cmap = colors.ListedColormap(['#ffdddd','#000000', '#ff9999', '#ff7777', '#ff5555', '#ff3333' ,'#ff1111'])
#'gray','green','pink','brown','blue','violet','yellow','orange', 
cmap = colors.ListedColormap(['#000000', 'gray','#ffdddd', '#ff9999', '#ff7777', '#ff5555', '#ff3333', '#ff1111', 'blue'])
boundaries = [np.power(10., 10), np.power(10., 11), np.power(10., 12), np.power(10., 13), np.power(10., 14), np.power(10., 15), np.power(10., 16), np.power(10.,17), np.power(10.,20.3), np.power(10.,25)]
norm = colors.BoundaryNorm(boundaries, cmap.N, clip=True)

pcmz0 = ax0.pcolormesh(xi, yi, z0, cmap=cmap, norm=norm)
ax0.set_xticks([x.min(), x.max()])
ax0.set_yticks([y.min(), y.max()])
ax0.set_xlim([x.min(), x.max()])
ax0.set_ylim([y.min(), y.max()])
ax0.set(aspect=1)
pcmz1 = ax1.pcolormesh(xi, yi, z1, cmap=cmap, norm=norm)
ax1.set_xticks([x.min(), x.max()])
ax1.set_yticks([y.min(), y.max()])
ax1.set_xlim([x.min(), x.max()])
ax1.set_ylim([y.min(), y.max()])
ax1.set(aspect=1)
pcmz2 = ax2.pcolormesh(xi, yi, z2, cmap=cmap, norm=norm)
ax2.set_xticks([x.min(), x.max()])
ax2.set_yticks([y.min(), y.max()])
ax2.set_xlim([x.min(), x.max()])
ax2.set_ylim([y.min(), y.max()])
ax2.set(aspect=1)

axparameters = plt.subplot(gs[1])
frame = plt.gca()
frame.axes.get_xaxis().set_ticks([])
frame.axes.get_yaxis().set_ticks([])
axparameters.set(aspect=1)
plt.text(0.1, 0.2, 'LAF=%.5f$\pm$%.5f \n\n\n LLS=%.5f$\pm$%.5f \n\n\n DLA=%.5f$\pm$%.5f' % (LAF, LAF_sigma, LLS, LLS_sigma, DLA, DLA_sigma), size=10)


axc = plt.subplot(gs[:,2])
plt.colorbar(pcmz1,  cax=axc,format=ticker.FuncFormatter(fmt))

ax0.set_ylabel("Y [kpc]")  
ax1.set_xlabel("X [kpc]")
ax1.set_ylabel("Z [kpc]") 
ax2.set_xlabel("Y [kpc]") 
fig.suptitle('HI Column Density Maps', size=20)
gs.tight_layout(fig, rect=[0, 0, 1, 0.9])
plt.savefig('m11q.png', dpi=900)
plt.show()

Which produces the following maps:

---

And here is the final version of the code which does everything I aimed for except few caveats mentioned inside the answer (by me) to the post:

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate
#import matplotlib.figure as fig
import matplotlib.cm as cm
import matplotlib.colors as colors
import matplotlib.colorbar as cb
import matplotlib.ticker as tk 
import pylab as pl
import matplotlib.ticker as ticker
from matplotlib import gridspec

x, y, covering_fraction_z, covering_fraction_x, covering_fraction_y = np.genfromtxt(r'm11q.txt', unpack=True)

def fmt(x, pos):
    return '10e{}'.format(np.log10(x))

nInterp = 3000
xi, yi = np.linspace(x.min(), x.max(), nInterp), np.linspace(y.min(), y.max(), nInterp)
xi, yi = np.meshgrid(xi, yi)

z0 = scipy.interpolate.griddata((x, y), covering_fraction_z, (xi, yi), method='nearest')
z1 = scipy.interpolate.griddata((x, y), covering_fraction_x, (xi, yi), method='nearest')
z2 = scipy.interpolate.griddata((x, y), covering_fraction_y, (xi, yi), method='nearest')


Pixel_Smoothing_Scale = 0.497
Halo_Pixels_No = (x.max()-x.min())*(y.max()-y.min())/ np.power(Pixel_Smoothing_Scale, 2)
DLA_Pixels_No_xy, DLA_Pixels_No_yz, DLA_Pixels_No_xz = 0, 0, 0
LLS_Pixels_No_xy, LLS_Pixels_No_yz, LLS_Pixels_No_xz = 0, 0, 0
LAF_Pixels_No_xy, LAF_Pixels_No_yz, LAF_Pixels_No_xz = 0, 0, 0

for i in np.arange(0, np.size(x)):
    if covering_fraction_z[i] >= np.power(10., 20.3):
        DLA_Pixels_No_xy+=1
    if ((covering_fraction_z[i] >= np.power(10., 17)) & (covering_fraction_z[i] < np.power(10., 20.3))):
        LLS_Pixels_No_xy+=1   
    if ((covering_fraction_z[i] >= np.power(10., 14)) & (covering_fraction_z[i] < np.power(10., 17))):
        LAF_Pixels_No_xy+=1    
    if covering_fraction_x[i] >= np.power(10., 20.3):
        DLA_Pixels_No_yz+=1
    if ((covering_fraction_x[i] >= np.power(10., 17)) & (covering_fraction_x[i] < np.power(10., 20.3))):
        LLS_Pixels_No_yz+=1 
    if ((covering_fraction_x[i] >= np.power(10., 14)) & (covering_fraction_x[i] < np.power(10., 17))):
        LAF_Pixels_No_yz+=1    
    if covering_fraction_y[i] >= np.power(10., 20.3):
        DLA_Pixels_No_xz+=1
    if ((covering_fraction_y[i] >= np.power(10., 17)) & (covering_fraction_y[i] < np.power(10., 20.3))):
        LLS_Pixels_No_xz+=1 
    if ((covering_fraction_y[i] >= np.power(10., 14)) & (covering_fraction_y[i] < np.power(10., 17))):
        LAF_Pixels_No_xz+=1

DLA_covering_fraction_XY, DLA_covering_fraction_YZ, DLA_covering_fraction_XZ = DLA_Pixels_No_xy/Halo_Pixels_No, DLA_Pixels_No_yz/Halo_Pixels_No, DLA_Pixels_No_xz/Halo_Pixels_No
DLA = np.sqrt(np.power(DLA_covering_fraction_XY, 2)+np.power(DLA_covering_fraction_YZ, 2)+np.power(DLA_covering_fraction_XZ, 2))/np.sqrt(3)
DLA_sigma = max([abs(DLA-DLA_covering_fraction_XY),abs(DLA-DLA_covering_fraction_YZ),abs(DLA-DLA_covering_fraction_XZ)])

LLS_covering_fraction_XY, LLS_covering_fraction_YZ, LLS_covering_fraction_XZ = LLS_Pixels_No_xy/Halo_Pixels_No, LLS_Pixels_No_yz/Halo_Pixels_No, LLS_Pixels_No_xz/Halo_Pixels_No
LLS = np.sqrt(np.power(LLS_covering_fraction_XY, 2)+np.power(LLS_covering_fraction_YZ, 2)+np.power(LLS_covering_fraction_XZ, 2))/np.sqrt(3)
LLS_sigma = max([abs(LLS-LLS_covering_fraction_XY),abs(LLS-LLS_covering_fraction_YZ),abs(LLS-LLS_covering_fraction_XZ)])

LAF_covering_fraction_XY, LAF_covering_fraction_YZ, LAF_covering_fraction_XZ = LAF_Pixels_No_xy/Halo_Pixels_No, LAF_Pixels_No_yz/Halo_Pixels_No, LAF_Pixels_No_xz/Halo_Pixels_No
LAF = np.sqrt(np.power(LAF_covering_fraction_XY, 2)+np.power(LAF_covering_fraction_YZ, 2)+np.power(LAF_covering_fraction_XZ, 2))/np.sqrt(3)
LAF_sigma = max([abs(LAF-LAF_covering_fraction_XY),abs(LAF-LAF_covering_fraction_YZ),abs(LAF-LAF_covering_fraction_XZ)])




fig = plt.figure(figsize=(8, 8))

gs = gridspec.GridSpec(2, 3, width_ratios=[5, 5, 1], height_ratios=[1,1,1]) 
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[3])
ax2 = plt.subplot(gs[4])

#cmap = colors.ListedColormap(['#ffdddd','#000000', '#ff9999', '#ff7777', '#ff5555', '#ff3333' ,'#ff1111'])
#'gray','green','pink','brown','blue','violet','yellow','orange', 
cmap = colors.ListedColormap(['#ffdddd', '#ff7777', '#ff3333', 'yellow', 'orange','green', 'red', 'blue'])
boundaries = [np.power(10., 4), np.power(10., 6), np.power(10., 8), np.power(10., 10), np.power(10., 12), np.power(10., 14), np.power(10.,17), np.power(10.,20.3), np.power(10.,25)]
norm = colors.BoundaryNorm(boundaries, cmap.N, clip=True)

pcmz0 = ax0.pcolormesh(xi, yi, z0, cmap=cmap, norm=norm)
ax0.set_xticks([x.min(), x.max()])
ax0.set_yticks([y.min(), y.max()])
ax0.set_xlim([x.min(), x.max()])
ax0.set_ylim([y.min(), y.max()])
ax0.set(aspect=1)
pcmz1 = ax1.pcolormesh(xi, yi, z1, cmap=cmap, norm=norm)
ax1.set_xticks([x.min(), x.max()])
ax1.set_yticks([y.min(), y.max()])
ax1.set_xlim([x.min(), x.max()])
ax1.set_ylim([y.min(), y.max()])
ax1.set(aspect=1)
pcmz2 = ax2.pcolormesh(xi, yi, z2, cmap=cmap, norm=norm)
ax2.set_xticks([x.min(), x.max()])
ax2.set_yticks([y.min(), y.max()])
ax2.set_xlim([x.min(), x.max()])
ax2.set_ylim([y.min(), y.max()])
ax2.set(aspect=1)

axparameters = plt.subplot(gs[1])
frame = plt.gca()
frame.axes.get_xaxis().set_ticks([])
frame.axes.get_yaxis().set_ticks([])
axparameters.set(aspect=1)
plt.text(0.1, 0.2, 'LAF=%.6f$\pm$%.6f \n\n\n LLS=%.6f$\pm$%.6f \n\n\n DLA=%.6f$\pm$%.6f' % (LAF, LAF_sigma, LLS, LLS_sigma, DLA, DLA_sigma), size=9)


axc = plt.subplot(gs[:,2])
plt.colorbar(pcmz1,  cax=axc,format=ticker.FuncFormatter(fmt))

ax0.set_ylabel("Y [kpc]")  
ax1.set_xlabel("X [kpc]")
ax1.set_ylabel("Z [kpc]") 
ax2.set_xlabel("Y [kpc]") 
fig.suptitle('HI Column Density Maps', size=20)
gs.tight_layout(fig, rect=[0, 0, 1, 0.9])
plt.savefig('m11q.png', dpi=900)
plt.show()

which produces the following maps:

Answer 1

(1) and (2) are related. You can define custom color scheme. For discrite arbitrary values it can be eg:

cmap = colors.ListedColormap(['red', '#000000','#444444', '#666666', '#ffffff', 'blue', 'orange'])
boundaries = [-1, -0.9, -0.6, -0.3, 0, 0.3, 0.6, 1]
norm = colors.BoundaryNorm(boundaries, cmap.N, clip=True)

This will map values following way:

• -1.0 x < -0.9 to red
• -0.9 < x < -0.6 to black (#000000)
• -0.6 < x < -0.3 to #444444
• and so on

Color i will be used for values between boundary i and i+1. Colors can be specified by names ('red', 'green'), HTML codes ('#ffaa44', '#441188') or RGB tuples ((0.2, 0.9, 0.45)).

Complete example:

# -*- coding: utf-8 -*-    
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate
#import matplotlib.figure as fig
import matplotlib.cm as cm
import matplotlib.colors as colors
import matplotlib.colorbar as cb
import matplotlib.ticker as tk 
import pylab as pl
import matplotlib.ticker as ticker
from matplotlib import gridspec

x, y, covering_fraction_z, covering_fraction_x, covering_fraction_y = np.genfromtxt(r'm11q.txt', unpack=True)

def fmt(x, pos):
    return '10e{}'.format(np.log10(x))

nInterp = 402
xi, yi = np.linspace(x.min(), x.max(), nInterp), np.linspace(y.min(), y.max(), nInterp)
xi, yi = np.meshgrid(xi, yi)

z0 = scipy.interpolate.griddata((x, y), covering_fraction_z, (xi, yi), method='nearest')
z1 = scipy.interpolate.griddata((x, y), covering_fraction_x, (xi, yi), method='nearest')
z2 = scipy.interpolate.griddata((x, y), covering_fraction_y, (xi, yi), method='nearest')

fig = plt.figure(figsize=(6, 6))

gs = gridspec.GridSpec(2, 3, width_ratios=[4, 4, 1], height_ratios=[1,1,1]) 
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[3])
ax2 = plt.subplot(gs[4])

cmap = colors.ListedColormap(['#ffdddd','#000000', '#ff9999', '#ff7777', '#ff5555', '#ff3333' ,'#ff1111'])
boundaries = [np.power(10., 14), np.power(10., 15.5), np.power(10.,17), np.power(10.,20.3), np.power(10.,25)]
norm = colors.BoundaryNorm(boundaries, cmap.N, clip=True)

pcmz0 = ax0.pcolormesh(xi, yi, z0, cmap=cmap, norm=norm)
ax0.set_xticks([x.min(), x.max()]) #arbitrary, also for Y axis ticks
ax0.set_xlim([x.min(), x.max()])
ax0.set_ylim([y.min(), y.max()])
pcmz1 = ax1.pcolormesh(xi, yi, z1, cmap=cmap, norm=norm)
ax1.set_xticks([x.min(), x.max()])
ax1.set_xlim([x.min(), x.max()])
ax1.set_ylim([y.min(), y.max()])
pcmz2 = ax2.pcolormesh(xi, yi, z2, cmap=cmap, norm=norm)
ax2.set_xticks([x.min(), x.max()])
ax2.set_xlim([x.min(), x.max()])
ax2.set_ylim([y.min(), y.max()])

axc = plt.subplot(gs[:,2])
plt.colorbar(pcmz1,  cax=axc,format=ticker.FuncFormatter(fmt))

ax0.set_ylabel("Y [kpc]")  
ax1.set_xlabel("X [kpc]")
ax1.set_ylabel("Z [kpc]") 
ax2.set_xlabel("Y [kpc]")   
plt.tight_layout()
plt.savefig('m11q.png', dpi=900)
plt.show()

For (3). There are many ways to do this. For example:

fraction = np.count_nonzero((np.power(10, 14) <= Z) & (Z < np.power(10, 15.5)) * 1.0 / Z.size

will give you fraction of values in range [10e14, 10e15.5)

(4) depends on method of creating subplots. With GridSpec it won't appear if you do not declare it explicitly.

Answer 2

Here is the answer to the updated version of the question. Originally I was not able to produce a color bar for my maps. Xevaquor helped me with that by the above answer. So, the answer was legitimately posted by Xevaquor and it is marked as "the answer" now. However, as I was taking my time to resolve some of the issues with covering fractions, finally I was able to understand the problem mathematically and then try to approach it numerically in a python code. My issue was that I did not realize that what I am after is in fact a ratio of the number of pixels. This is why, I had to count the number of pixels that lie inside a range of column densities and then divide that by the total number of pixels making each projection of the map. Since the smoothing scale is known to be $0.497 kpc$, we know that total number of pixels will be the physical area of the map in $kpc^{2}$ divided by the area of each pixel namely $0.497^{2}$. So, the following snippet of code can be added to the original and/or Xevaquor's code to account for the covering fractions correctly:

Pixel_Smoothing_Scale = 0.497
Halo_Pixels_No = (x.max()-x.min())*(y.max()-y.min())/ np.power(Pixel_Smoothing_Scale, 2)
DLA_Pixels_No_xy, DLA_Pixels_No_yz, DLA_Pixels_No_xz = 0, 0, 0
LLS_Pixels_No_xy, LLS_Pixels_No_yz, LLS_Pixels_No_xz = 0, 0, 0

for i in np.arange(0, np.size(x)):
    if covering_fraction_z[i] >= np.power(10., 20.3):
        DLA_Pixels_No_xy+=1
    if ((covering_fraction_z[i] >= np.power(10., 17)) & (covering_fraction_z[i] < np.power(10., 20.3))):
        LLS_Pixels_No_xy+=1   
    if covering_fraction_x[i] >= np.power(10., 20.3):
        DLA_Pixels_No_yz+=1
    if ((covering_fraction_x[i] >= np.power(10., 17)) & (covering_fraction_x[i] < np.power(10., 20.3))):
        LLS_Pixels_No_yz+=1 
    if covering_fraction_y[i] >= np.power(10., 20.3):
        DLA_Pixels_No_xz+=1
    if ((covering_fraction_y[i] >= np.power(10., 17)) & (covering_fraction_y[i] < np.power(10., 20.3))):
        LLS_Pixels_No_xz+=1 

DLA_covering_fraction_XY, DLA_covering_fraction_YZ, DLA_covering_fraction_XZ = DLA_Pixels_No_xy/Halo_Pixels_No, DLA_Pixels_No_yz/Halo_Pixels_No, DLA_Pixels_No_xz/Halo_Pixels_No
DLA = np.sqrt(np.power(DLA_covering_fraction_XY, 2)+np.power(DLA_covering_fraction_YZ, 2)+np.power(DLA_covering_fraction_XZ, 2))/np.sqrt(3)
DLA_sigma = max([abs(DLA-DLA_covering_fraction_XY),abs(DLA-DLA_covering_fraction_YZ),abs(DLA-DLA_covering_fraction_XZ)])

LLS_covering_fraction_XY, LLS_covering_fraction_YZ, LLS_covering_fraction_XZ = LLS_Pixels_No_xy/Halo_Pixels_No, LLS_Pixels_No_yz/Halo_Pixels_No, LLS_Pixels_No_xz/Halo_Pixels_No
LLS = np.sqrt(np.power(LLS_covering_fraction_XY, 2)+np.power(LLS_covering_fraction_YZ, 2)+np.power(LLS_covering_fraction_XZ, 2))/np.sqrt(3)
LLS_sigma = max([abs(LLS-LLS_covering_fraction_XY),abs(LLS-LLS_covering_fraction_YZ),abs(LLS-LLS_covering_fraction_XZ)])

print(Pixel_Smoothing_Scale)
print(LLS_covering_fraction_XY, LLS_covering_fraction_YZ, LLS_covering_fraction_XZ)
print(DLA_covering_fraction_XY, DLA_covering_fraction_YZ, DLA_covering_fraction_XZ)

---

With this being said, the only issue that is unresolved is the slowness of the code for some files that are huge and that increasing interp value cannot be continued indefinitely towards larger values in order to get smoother maps. So, I don't know how to approach this problem in order to make smoother maps (colors be filled homogeneously and no single region on the map be devoid of any color) without making the code increasingly slow. Any help towards this issue is appreciated.