This is the first part - building a polygon grid from your center coordinates, and outputting it with a few basic attributes to shapefile. The code I've copied below is a script based on the QGIS ftools VectorGrid function, which basically just lets you build a polygon grid and write to shapefile without the QGIS environment.
For the second part - counting the data points within each grid - you can either add some lines within the Mk_polygrid function which find datapoints within each polygon as it is being built and written, or you can take the output shapefile and add fields counting your points in another step.
the well-known-text I could find for RT90 2.5 gon W (link: http://epsg.io/2400 ) looks like this, but notes that it is deprecated:
rt90 = 'PROJCS["RT90 2.5 gon W (deprecated)",GEOGCS["RT90",DATUM["Rikets_koordinatsystem_1990",SPHEROID["Bessel1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[414.1,41.3,603.1,-0.855,2.141,-7.023,0],AUTHORITY["EPSG","6124"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4124"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15.80827777777778],PARAMETER["scale_factor",1],PARAMETER["false_easting",1500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","2400"]]'
import ogr,sys,math, osr
params = {}
### DEFINE PROJECTION
# Replace the projection string with your own (RT90 2.5 GW)
# rt90 = 'ogc well-known-text projection string'
params['prj'] = rt90
### DEFINE OUTFILE AND PROJECTION
params['outf'] = 'yourfile'
# no extension; script will make .shp, .shx, .dbf,and .prj files
### DEFINE YOUR GRID SPACING
params['dx'] = 100 #meters
params['dy'] = 100 #meters
# assuming x and y coords are spaced as multiples of 100:
params['minx'] = min(your_x_data)
params['maxx'] = max(your_x_data)
params['miny'] = min(your_y_data)
params['maxy'] = max(your_y_data)
### DEFINE YOUR GRID BUILDING FUNCTION
def Mk_polygrid(minx,maxx,miny,maxy,dx,dy,prj,outf):
'''QGIS fTools approach (VectorGrid tool);
modified to replace QGIS stuff with OGR
and my little extent grabber functions'''
driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(outf+'.shp')
layer = ds.CreateLayer('',None,ogr.wkbPolygon)
layer.CreateField(ogr.FieldDefn('id',ogr.OFTInteger))
layer.CreateField(ogr.FieldDefn('ctr_x',ogr.OFTReal))
layer.CreateField(ogr.FieldDefn('ctr_y',ogr.OFTReal))
defn = layer.GetLayerDefn()
idVar = 0
progVar = 0
count = 0
count_max = (maxy-miny) / dy
count_update = count_max * 0.05 # print progress every 5%!
y = maxy
while int(y) > int(miny):
x = minx
while int(x) < int(maxx):
perim = ogr.Geometry(ogr.wkbLinearRing)
perim.AddPoint(x - 0.5*dx, y - 0.5*dy)
perim.AddPoint(x + 0.5*dx, y - 0.5*dy)
perim.AddPoint(x + 0.5*dx, y + 0,5*dy)
perim.AddPoint(x - 0.5*dx, y + 0.5*dy)
perim.AddPoint(x - 0.5*dx, y - 0.5*dy)
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(perim)
feat = ogr.Feature(defn)
feat.SetField('id',idVar)
feat.SetField('ctr_x',x)
feat.SetField('ctr_y',y)
feat.SetGeometry(polygon)
layer.CreateFeature(feat)
feat = geom = None
idVar += 1
x += dx
y -= dy
count += 1
if int( math.fmod( count, count_update ) ) == 0:
prog = int( count / count_max * 100 )
report = '%s%% . . ' % prog
sys.stdout.write( report )
sys.stdout.flush()
# Save and close everything
ds = layer = feat = perim = polygon = None
sys.stdout.write("\n")
### DEFINE YOUR PRJ BUILDING FUNCTION
def Mk_proj(prj,outf):
with open("%s.prj" % outf, "w") as proj:
proj.write(prj)
### SCRIPT ACTIONS
if __name__ == '__main__':
Mk_polygrid(**params)
Mk_proj(params['prj'],params['outf'])