pyshp not working properly when creating polygons (offset + missing polygons)

Question

I have a esri shapefile thats projected in NSIDC EASE-Grid Global that has roughly 50k polygons. I wrote the python script below to create a grid that covers the span of each polygons bbox (note I don't care if there are rectangles that don't interesect a polygon because I am going to use QGIS to remove any that dont, my computer just doesn't have enough memory/disc space to create a grid that covers the extent of the whole world).

However, my code seems to have an offset i.e all the rectangles are off by a constant ammount from where they should be. Furthermore, the script just doesn't create a grid of rectangles for some polygons in the shapefile??? I have checked the validity of all the polygons and they all check out. I'm not sure whats going wrong... can someone help me figure this out?

import shapefile
import math
import time
sf = shapefile.Reader("2500_mines_NSDIC_EASE/2500_mines_NSDIC_EASE.shp")

shapes = sf.shapes()
w = shapefile.Writer('grid/grid', shapeType=shapefile.POLYGON)
w.field('name', 'C')
count = 1
start = time.time()
def create_poly(bottom_x, bottom_y, width, height):
    global count
    w.poly([[[bottom_x, bottom_y],[bottom_x+width, bottom_y], [bottom_x+width, bottom_y+height], [bottom_x, bottom_y+height], [bottom_x, bottom_y]]])
    w.record('rect' + str(count))
    count += 1
tile_size = 400#pixels
resolution = 10 #m/pixel
rect_dim = tile_size * resolution
num_processed = 0
for shape in shapes:
    bbox = shape.bbox
    minx = bbox[0]
    miny = bbox[1]
    maxx = bbox[2]
    maxy = bbox[3]
    assert(minx < maxx)
    assert(miny < maxy)

    minx_tile = math.floor(minx/rect_dim) -1
    miny_tile = math.floor(miny/rect_dim) -1
    maxx_tile = math.ceil(maxx/rect_dim) + 1
    maxy_tile = math.ceil(maxy/rect_dim) + 1

    for x in range(minx_tile, maxx_tile+1):##inclusive, exclusive
        for y in range(miny_tile, maxy_tile+1):
            create_poly(x*rect_dim, y*rect_dim, rect_dim, rect_dim)
    num_processed += 1
    if (num_processed % 1000 == 0):
        print("Processed " + str(num_processed) + " mines. " + str(count) + " polygons " + str(time.time()-start))
        

w.close()

Heres an example of the offset (the other rectangles for other polygons are offset by the same ammount)(viewed in QGIS)

For some reason there appears to be no overlap in Northern Africa?