1

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?

  • Do you want the result output as a shapefile or another GIS format, or do you want to keep it in memory then do stuff with the group of polygons? – WileyB May 20 '14 at 9:48
  • @WileyB I want it as a shapefile. I will later on 'crop' the grid using another shapefile with administrative regions to keep the squares inside the city I'm mapping. – Kenji May 20 '14 at 10:14
  • OK. Do you know if your coordinate system is a well known / common spatial reference system (like UTM) or listed in the EPSG registry? – WileyB May 20 '14 at 10:22
  • @WileyB It is the RT90 2.5 GW – Kenji May 20 '14 at 10:26
3

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'])
  • That is much more help than I had expected. Thank you so much for this and for reminding me that RT90 is deprecated. I should have expected so as my data is old. – Kenji May 20 '14 at 11:13
  • 1
    RT90 2.5 gon V seems to be identical in the epsg registry, and not deprecated: epsg.io/3021 – WileyB May 20 '14 at 11:14
  • 1
    Bizarre coincidence; I was working on these functions earlier this week, and even partly with the same projection! – WileyB May 20 '14 at 11:16
  • I guess from the moose in your profile that we are in the same part of the world. Thank you again for the help. – Kenji May 20 '14 at 11:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.