Mathematically producing sphere-shaped hexagonal grid

The shape you have is one of so called “Goldberg polyhedra”, is also a geodesic polyhedra.

The (rather elegant) algorithm to generate this (and many many more) can be succinctly encoded in something called a Conway Polyhedron Notation.

The construction is easy to follow step by step, you can click the images below to get a live preview.

1. The polyhedron you are looking for can be generated from an icosahedron — Initialise a mesh with an icosahedron.
   

2. We apply a “Truncate” operation (Conway notation t) to the mesh (the sperical mapping of this one is a football).
   

3. We apply the “Dual” operator (Conway notation d).
   

4. We apply a “Truncate” operation again. At this point the recipe is tdtI (read from right!). You can already see where this is going.
   

5. Apply steps 3 & 4 repeatedly until you are satisfied.

For example below is the mesh for dtdtdtdtI.

This is quite easy to implement. I would suggest using a datastructure that makes it easy to traverse the neighbourhood give a vertex, edge etc. such as winged-edge or half-edge datastructures for your mesh. You only need to implement truncate and dual operators for the shape you are looking for.