Convert lines density into meshgrid (Python to GMT)

This might be a very basic question, I have plot the 2d ray path but Can we do the following work in GMT6? Any suggestion would be appreciated

I have multiple lines within region where x= 90.3 to 91 and y= 20.3 to 20.7.
I want to divide the map into smaller grid sizes (0.001) to convert like points and plot the point density (line) in each grid using colormesh. I have done some tests in Python using the below dataset but it seems like it also shows the density where no lines are crossing. I also wanna overlap some feature in GMT.

import numpy as np
import matplotlib.pyplot as plt

# map range and grid size
x_min, x_max = 90.4, 90.9
y_min, y_max = 20.3, 20.7
grid_size = 0.001

# grid of x and y coordinates
x = np.arange(x_min, x_max, grid_size)
y = np.arange(y_min, y_max, grid_size)
X, Y = np.meshgrid(x, y)

density_grid = np.zeros_like(X)

# demo lines (start and end points)
lines = [
[(90.42, 20.32), (90.45, 20.6)],
[(90.48, 20.31), (90.69, 20.32)],
[(90.47, 20.5), (90.89, 20.42)],
[(90.41, 20.33), (90.45, 20.5)],
[(90.43, 20.34), (90.47, 20.48)],
]

# Iterate through each line and update density grid
for line in lines:
(x1, y1), (x2, y2) = line
indices = np.where((x1 <= X) & (X <= x2) & (y1 <= Y) & (Y <= y2))
density_grid[indices] += 1

# Plotting the map, density grid, and lines
fig, ax = plt.subplots()
ax.scatter(X, Y, color='black', marker='.')
mesh = ax.pcolormesh(X, Y, density_grid, cmap='hot', shading='auto')
for line in lines:
(x1, y1), (x2, y2) = line
ax.plot([x1, x2], [y1, y2], color='blue')

plt.colorbar(mesh, ax=ax)
plt.show()

Please, use the PyGMT category to post PyGMT questions.

1 Like