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.