Access color of a point given the 3D position (on the surface of a polygon)
Some indications for a way to do that (maybe there are other possibilities... I'd be interested if some other answer could offer something more simple).

From the comments below the question, I assume the following entry parameters:


A mesh
A ray cast point 'V' on the mesh, so you also have the corresponding polygon index 'P' (considering you are using BVHTree from the mathutils module)
A UV map
A texture


First thing: at some point, you'll need to match V position on P to the V' position on the UV map part where P correspond to.
That can be done using barycentric_transform (function which can be found in mathutils.geometry), but a triangle is needed (as we need all to be convex).
This means that you need to work on a triangulated mesh (you could also triangulate on the fly, but I don' know how to do it).

So the overall process can be the following:


Get your mesh
Triangulate it
Ray cast and obtain V and P (a triangle)
Find the UV map part P' corresponding to P
Calculate V' using barycentric_transform
Remap V' on the texture


Some indications for all that using bpy and bmesh modules


Get your mesh


Snippet

import bpy
obj = bpy.data.object['your mesh name']



Triangulate it


Snippet

import bpy
import bmesh
def TriangulateMesh( obj ):
    bm = bmesh.new()
    bm.from_mesh( obj.data )
    bmesh.ops.triangulate( bm, faces=bm.faces[:] )
    bm.to_mesh( obj.data )
    bm.free()


Eventually, work on a copy:

def CopyMesh( obj ):
    copy = obj.copy()
    copy.data = obj.data.copy()
    bpy.context.scene.objects.link( copy )
    return copy



Ray cast


Snippet

from mathutils.bvhtree import BVHTree
... some code and at the end:
location, normal, index, distance = bvh.ray_cast( ... your parameters )


index here is P.


Find the UV map part corresponding to P


Get 3D vertices indices

verticesIndices = obj.data.polygons[P].vertices


P vertices coordinates in 3D space

p1, p2, p3 = [obj.data.vertices[verticesIndices[i]].co for i in range(3)]


P UV map vertices indices are found in:

uvMapIndices = obj.data.polygons[P].loop_indices


The UV map

uvMap = obj.data.uv_layers['your uv map name']


The coordinates of P vertices in the UV map space:

uv1, uv2, uv3 = [uvMap.data[uvMapIndices[i]].uv for i in range(3)]



Calculate V' using barycentric_transform


V' is given by (make uv1, 2 and 3 as 3D vectors before):

V' = barycentric_transform( V, p1, p2, p3, uv1, uv2, uv3 )



Remap V' on the texture


Get the texture, the pixels and texture dimensions:

image = bpy.data.images['your texture name']
pixels = list( image.pixels ) #Faster than accessible image.pixels[x] each time
width = image.size[0]
height = image.size[1]


From V', get x and y. These values as give as if the image size was (1, 1). So you need to rescale and modulo x and y, then reach the corresponding image color using 'pixels' (which is a linear byte array of individual color components RGBA).