python - Creating counting grid from set of points

Question

I have a dataset of longitudes/latitudes as follows:

id,spp,lon,lat
1a,sp1,1,9
1b,sp1,3,11
1c,sp1,6,12
2a,sp2,1,9
2b,sp2,1,10
2c,sp2,3,10
2d,sp2,4,11
2e,sp2,5,12
2f,sp2,6,12
3a,sp3,4,13
3b,sp3,5,11
3c,sp3,8,8
4a,sp4,4,12
4b,sp4,6,11
4c,sp4,7,8
5a,sp5,8,8
5b,sp5,7,6
5c,sp5,8,2
6a,sp6,8,8
6b,sp6,7,5
6c,sp6,8,3

From such data, I want to generate a grid like this:

0 0 0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 1 0 0 0 0 0 0 
0 0 0 0 0 0 1 1 2 0 0 0 0 
0 0 0 0 0 1 1 1 1 0 0 0 0 
0 0 0 1 0 1 0 0 0 0 0 0 0 
0 0 0 2 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 1 3 0 0 
0 0 0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 1 0 0 
0 0 0 0 0 0 0 0 0 0 1 0 0 
0 0 0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 0 0

which gives the number of data records in each cell of the grid, using variable "spp" as a categorical (grouping) factor.

From this grid, I then want to create a heat map, superimposed on a geographical map, so that I end up with something like the figure below.

I (finally) managed to write some code which does what I want.

Here it is:

import csv
import numpy as np
from mpl_toolkits.basemap import Basemap
from matplotlib import cm as cmap
import matplotlib.pyplot as plt

import warnings
warnings.filterwarnings("ignore")

#read input data
with open('testdata.csv', 'r') as csvfile:
    reader = csv.reader(csvfile, delimiter=',')
    headers = reader.next()
    input_data = list(reader)

#grid dimensions with one-degree resolution
lat_inf, lon_inf = 0, 0
lat_sup, lon_sup = 90, 360
resolution = 1

latitude, longitude = [], []
latitude = range(lat_inf, lat_sup, resolution)
longitude = range(lon_inf, lon_sup, resolution)

#create output grid
output_grid = []
for i in latitude:
   output_grid.append([])
   for j in longitude:
       output_grid[i].append(0)

#traverse the input_data evaluating the lat, lon coordinates
#summing +1 in the output_grid[latitude][longitude].
for row in input_data:
   lat = int(row[2]) 
   lon = int(row[3]) 
   #sp = row[1]

   #check its indexes
   i_lat = latitude.index(lat)
   i_lon = longitude.index(lon)

   #increase counter
   output_grid[i_lat][i_lon] += 1

output_grid = np.array(output_grid, np.int16)

#create map
m = Basemap()
m.drawcoastlines(linewidth=0.25)

#display image
im = m.imshow(output_grid.transpose(), cmap='summer', origin='lower', aspect='auto', interpolation='none')
m.colorbar(im)
plt.show()

It (mostly) works, but te problem is that the grid image is not being displayed correctly: it appears too small, in the lower left corner of the map, as shown in the figure below).

Also, is there a way to change the backgound color of the image grid, other than fiddling with Matplotlib colormaps?

Any hints, ideas, suggestions will be much appreciated.

Answer 1

Drawing an image on a basemap might be possible, but you would need to fiddle with the coordinates and the extent to get it at the right position.

Better use a pcolormesh which supports non-equally spaced grids and hence is well suited for different projections. What you need is then not only your data, but also a grid. This grid defines the location of each point of the data. Before using the grid, you need to transform it to the basemap coordinate system, so if lons,lats = np.meshgrid(..., ...) is the grid, and m = Basemap(..), the transformed grid is

X, Y = m(lons,lats)

This can then be supplied to pcolormesh

m.pcolormesh(X,Y,data)

A complete example, where we have a data array of shape (60x40) which we want to reside in the longitude range between -10 and 50 degrees and in the latitude range between -20 and 20 degrees:

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap

data = np.random.rand(40,60)

lons,lats = np.meshgrid(np.linspace(-10,50,data.shape[1]+1),
                        np.linspace(-20,20,data.shape[0]+1))

fig, ax = plt.subplots()
m = Basemap(projection='merc', ax=ax, lat_0=0.,lon_0=0.,
            llcrnrlon=-179.,llcrnrlat=-80.,urcrnrlon=170.,urcrnrlat=80.)
m.drawcoastlines()
X, Y = m(lons,lats)
pc = m.pcolormesh(X,Y,data,cmap='RdBu_r')

m.drawparallels(np.arange(-60,61,20),labels=[1,1,0,1])
m.drawmeridians([-90,-10,50,90],labels=[1,1,0,1])
plt.show()

Answer 2

The doc for Basemap says

Calling a Basemap class instance with the arguments lon, lat will convert lon/lat (in degrees) to x/y map projection coordinates (in meters). The inverse transformation is done if the optional keyword inverse is set to True.

So, you need to covert your data to degrees of lat and lon, pass those to Basemap, and use what is returned to you to plot your locations.

The color you're getting is the default color for "cool" (see the color bar to the right of your plot). I'm not sure how to change it, but I think there should be arguments you can pass to MPL routines to tell it which color to use for "cool" and which to use for "warm".