Creating grid polygons from coordinates using R or Python

I have locational data with north-south and east-west coordinates, with a resolution of 100 meters squares. Each data point's location is somewhere within the 100 meters square which centre is given by the coordinates. The values of the north-south coordinates show the distance from the equator, and the values of east-west coordinates show the distance from the meridian 1,500,000 meters east of the Greenwich meridian. The coordinate system is the RT90 2.5 GW.

I want to create a square grid polygon from these data, with the number of data points and the centre's coordinates as an attribute of each square in the grid. I need the centres of the squares in the grid to match the coordinates in the dataset.

I would like to have the result as a shapefile, as I will 'crop' the grid with another shapefile for administrative regions of the city I'm mapping. Is there a way of doing it using either R or Python?

Answer

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'])