For high quality work, the best approach is to fit a [generalized linear spatial model](http://cran.r-project.org/web/packages/geoRglm/index.html) for a multinomial process.  That's a procedure that requires considerable expertise and effort.  As a substitute, you might consider **expanding each sample point into its polygon of influence** (aka Thiessen polygon, Voronoi polygon, or Dirichlet cell).  Limiting the expansion to land areas is a good idea; this can be done with a mask grid.

**To illustrate,** consider this much smaller dataset (of 14,136 points) representing 12 lithologic classes as distinguished by color:

![Samples][1]

I accomplished the expansion by converting these points into a grid (around 800 rows and 1000 columns) and computing their [Euclidean allocation](http://edndoc.esri.com/arcobjects/9.2/net/shared/geoprocessing/spatial_analyst_tools/euclidean_allocation.htm), using a mask that limited the calculation to non-glaciated land.  (The color scheme in the next two figures differs from that of the preceding one.)

![Euclidean allocation][2]

For comparison, here is a [detailed lithologic map](http://projects.gtk.fi/export/sites/projects/barents/images/soil2.gif) of the same area drawn to the same scale with the same symbolization:

![Original map][3]

With a truly large dataset or a convoluted study area, it may be expeditious to tile the region and perform this procedure separately on each tile, mosaicing the results into one output raster if desired.  For this to work, the tiles need to overlap slightly to avoid edge effects (and then should be uniformly trimmed before mosaicing).

**The principal reasons for going to a raster representation** are (1) it's quick and easy to compute and (2) accurate vector-based solutions will be hard to come by.  If you try buffers, convex hulls, concave hulls, or whatever, you will find that they all mutually intersect and they still leave gaps: in other words, they won't produce a topologically consistent partition of the space into distinct lithological domains.

One vector-based method that *will* work is to compute a constrained Voronoi tessellation of the points (good methods take O(n*log(n)) time for n points), spatially merge the Voronoi cells according to the lithological attributes of their associated points, and then separate the resulting multi-polygons into their connected components (if you wish).  However, if all you need is vector *output*, it's easier to regiongroup the raster result and convert that to vector format.


  [1]: https://i.sstatic.net/oPu6d.png
  [2]: https://i.sstatic.net/wCYMD.png
  [3]: https://i.sstatic.net/g5MVo.png