Splitting a Shapefile based off coordinates in Python using osgeo (ogr) or pyshp

Question

I'm having trouble with splitting a shapefile of all the regions in Mongolia into a smaller shapefile based off some input coordinates.

The big shapefile I'm trying to split can be found here: https://marine.rutgers.edu/~cfree/gis-data/mongolia-gis-data/ under Administrative areas - soum boundaries.

I currently have input coordinates of upper left corner:

  uly, ulx = 47.82259512683848,108.3279267594774 

and bottom right corner:

  lry, lrx = 46.85958308337844,109.81664846692111

I've tried importing the shapefile with pyshp and the ogr, and have been having the most success with pyshp, using:

import shapefile
sf = shapefile.Reader('file path')
sf.bbox()

The results of sf.bbox() return [-325212.4151875854, 4620806.467562388, 2096670.9522781507, 5777542.4030963285], and I've been trying to find out what projection this is in, I haven't found a clear way to do this with pyshp.

Once I've worked out the current projection of the shapefile and reprojected it into lat/long coordinates, I then need to split it to get rid of everything outside the lat/long coordinates described above.

I've tried to split it so far using the following code, but this produced a shapefile containing nothing.

w = shapefile.Writer(sf.shapeType)
for i in sf.iterShapeRecords():
    geom = i.shape
    rec = i.record
    sxmin, symin, sxmax, symax = geom.bbox
    if sxmin < ulx: continue
    elif sxmax > lrx: continue
    elif symin < uly: continue
    elif symax > lry: continue
    w._shapes.append(geom)
    w.records.append(rec)
w.save('Local File')

I can try to provide more details if required.

Answer 1

Right, so after calming down and coming back to this I've managed to solve it, using the code below and opening the .prj in text editor to find the UTM projection of the shapefile.

import shapefile
sf = shapefile.Reader('file path')
w = shapefile.Writer(sf.shapeType)
w.fields = list(sf.fields)

source = osr.SpatialReference()
source.ImportFromEPSG(32647)

target = osr.SpatialReference()
target.ImportFromEPSG(4326)
transform = osr.CoordinateTransformation(source,target)

for layer in sf.iterShapeRecords():
    rec = layer.record
    xmin, ymax, xmax, ymin = layer.shape.bbox
    xmin, ymin, zmin = transform.TransformPoint(xmin,ymin)
    xmax, ymax, zmax = transform.TransformPoint(xmax,ymax)
    if (longs[0] < xmin < longs[1]) & (lats[0] > ymin > lats[1]):
        w._shapes.append(layer.shape)
        w.records.append(rec)
    if (longs[0] < xmax < longs[1]) & (lats[0] > ymax > lats[1]):
        w._shapes.append(layer.shape)
        w.records.append(rec)
w.save('clipped')