Create an UTM/MGRS graticule/fishnet in Python

Question

Following Creating square grid polygon shapefile with Python?, I'd like to create an UTM or MGRS graticule in python. In the mentioned question the given answer gives me a way to generate a grid but I don't know if the SRS or the projection is wrong, but the grid is not aligned with the graticule available on NGA.

Here is what I get :

Red grid is the grid from NGA, it aligns well in regard with UTM zone 31 and distance is good too. Blue grid is the one generated by the script I got from Creating square grid polygon shapefile with Python?, but modified to take geographic bounds, and intended utm zone. There is a slight offset (better than previously)

import os, sys
import ogr
import osr
from math import ceil


def main(outputGridfn,lonmin,lonmax,latmin,latmax, utmZone, gridHeight,gridWidth):

   # create the spatial reference, WGS84
    targetSRS = osr.SpatialReference()


    targetSRS.SetProjCS( "UTM "+ utmZone+ " (WGS84) in northern hemisphere." );
    targetSRS.SetWellKnownGeogCS( "WGS84" );
    targetSRS.SetUTM( int(utmZone), True );

    sourceSRS  = osr.SpatialReference()
    sourceSRS.ImportFromEPSG(4326)

    poCT =osr.CoordinateTransformation ( sourceSRS, targetSRS )

    # convert sys.argv to float
    lonmin = float(lonmin)
    lonmax = float(lonmax)
    latmin = float(latmin)
    latmax = float(latmax)

    xmin,ymin,_ = poCT.TransformPoint(lonmin,latmin)
    print (xmin, ymin)
    xmax,ymax,_ = poCT.TransformPoint(lonmax,latmax)

    gridWidth = float(gridWidth)
    gridHeight = float(gridHeight)

    # get rows
    rows = ceil((ymax-ymin)/gridHeight)
    # get columns
    cols = ceil((xmax-xmin)/gridWidth)

    # start grid cell envelope
    ringXleftOrigin = xmin
    ringXrightOrigin = xmin + gridWidth
    ringYtopOrigin = ymax
    ringYbottomOrigin = ymax-gridHeight

    # create output file
    outDriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputGridfn):
        os.remove(outputGridfn)
    outDataSource = outDriver.CreateDataSource(outputGridfn)


    outLayer = outDataSource.CreateLayer(outputGridfn,targetSRS,geom_type=ogr.wkbPolygon )

    featureDefn = outLayer.GetLayerDefn()

    # create grid cells
    countcols = 0
    while countcols < cols:
        countcols += 1

        # reset envelope for rows
        ringYtop = ringYtopOrigin
        ringYbottom =ringYbottomOrigin
        countrows = 0

        while countrows < rows:
            countrows += 1
            ring = ogr.Geometry(ogr.wkbLinearRing)
            ring.AddPoint(ringXleftOrigin, ringYtop)
            ring.AddPoint(ringXrightOrigin, ringYtop)
            ring.AddPoint(ringXrightOrigin, ringYbottom)
            ring.AddPoint(ringXleftOrigin, ringYbottom)
            ring.AddPoint(ringXleftOrigin, ringYtop)
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)

            # add new geom to layer
            outFeature = ogr.Feature(featureDefn)
            outFeature.SetGeometry(poly)
            outLayer.CreateFeature(outFeature)
            outFeature.Destroy

            # new envelope for next poly
            ringYtop = ringYtop - gridHeight
            ringYbottom = ringYbottom - gridHeight

        # new envelope for next poly
        ringXleftOrigin = ringXleftOrigin + gridWidth
        ringXrightOrigin = ringXrightOrigin + gridWidth

    # Close DataSources
    outDataSource.Destroy()


if __name__ == "__main__":

    #
    # example run : $ python grid.py <full-path><output-shapefile-name>.shp xmin xmax ymin ymax gridHeight gridWidth
    #

    if len( sys.argv ) != 9:
        print ("[ ERROR ] you must supply eight arguments: output-shapefile-name.shp lonmin lonmax latmin latmax utmZone gridHeight gridWidth")
        sys.exit( 1 )

    main( sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7],sys.argv[8] )

I'm quite confused at this time. Should I use directly UTM or MGRS coordinates?

I've seen in Generating shapefile containing UTM grid (graticule)? the option of using Saga GIS, but if possible I'd rather use code to be able to script it.

Edit: improved the computation