You can use the GDAL Python bindings. Examples on how to use it can you find here.
For example you create points with lat/lon like this
from osgeo import ogr # first import the library
point1 = ogr.Geometry(ogr.wkbPoint)
point1.AddPoint(13.381348,52.536273) # Berlin
point2 = ogr.Geometry(ogr.wkbPoint)
point2.AddPoint(11.557617,48.136767) # Munich
Create a transformation from EPSG:4326 (lat/long) to EPSG:3035 (Projected coordinate system for Europe in meters)
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(4326)
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(3035)
coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
Transform our points
point1.Transform(coordTransform)
>> POINT (4550348.379724434576929 3275002.535206703934819 0)
point2.Transform(coordTransform)
>>POINT (4436977.705661337822676 2781579.793173507787287 0)
And get the distance like this
point1.Distance(point2) # Distance in meter from Munich to Berlin
>>> 506279.480221 # roughly 506 km
Or you can create a polygon
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(-23.378906,68.974164) # North West Corner of Europe
ring.AddPoint(-23.378906,34.307144) # South West Corner of Europe
ring.AddPoint(31.464844,34.307144) # South East Corner of Europe
ring.AddPoint(31.464844,68.974164) # North East Corner of Europe
ring.AddPoint(-23.378906,68.974164) # North West Corner of Europe (to close to polygon)
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(ring) # rough bounding box around Europe
polygon.Transform(coordTransform) # transform it
And check if the polygon contains the point
polygon.Contains(point1) # Does Europe contain Berlin?
>>> True # It does ;)
And for the rest you can write functions to do what you want.