How to get the administrative region (admin_level=4) a place belongs to?

You will save yourself a lot of trouble if you can afford to run a local "Nominatim" installation and then simply use Nominatim's reverse geocoding ability. Alternatively, try "Photon" for the same purpose; there are ready-to-go Photon data dumps available for download, obviating the need for a database and an import process.

Otherwise, the likely easiest path is to pre-filter a planet file with the command-line osmium utility (only keep stuff tagged boundary=administrative and with an admin_level of interest to you), then load the entirety of that file into Pyosmium, construct geometries, use the Pyosmium/Shapely interoperability to stuff the geometries into a Shapely R-Tree. Once that is done, you can efficiently query which administrative bounding boxes a point is in; if this yields more than one result (eg half the points in New Mexico are also in the Texas bounding box), you need to do proper (and more costly) point-in-polygon checks for the resulting polygons but that is also easy with Shapely.

A third option is to do essentially what I wrote above but employ a PostGIS database to do it, i.e. load admin boundaries into PostGIS with osm2pgsql, then run SQL queries to find which polygons a point is in.