Dividing Geographical region into equal sized grid and retrieving indexing position

I have GPS coordinates of human mobility in Beijing Area. I want to divide geographical space in to rectangular grid, for example 2 sq. km (delta), and access the indexing position of any point inside the grid. The cells need not be exactly equal size, approximation can work in my case.

I have geographical region with the following coordinates (latitude, longitudes) of bounding box.

Bottom Left  (x1,y1) = 39.77750000, 116.17944444  
Top Left     (x1,y2) = 40.04722222, 116.17944444  
Bottom Right (x2,y1) = 39.77750000, 116.58888889  
Top Right    (x2,y2) = 40.04722222, 116.58888889  

This is rectangular region of 30 km x 34 kms.
The solution in my mind is to take delta as 2km, increment latitude and longitude value by delta until upper bound reaches.
To access the index position of GPS point p, Let BL be the Bottom Left point of rectangular region

Row =    Distance [(p.lat,BL.long), (BL.lat, BL.long)] / delta 
Column = Distance [(BL.lat,p.long), (BL.lat, BL.long)] / delta 

Is there any easier way or supporting library to solve this problem? preferably in combination of row and columns (x,y), so that I can measure the closeness of grid cells by finding distance between two grid cells in Cartesian coordinates system. An example image and input dataset may give you clear idea of the description.

The input dataset in given in the link https://drive.google.com/file/d/1JjvS7igTmrtLA4E5Rs5D6tsdAXqzpYqX/view