1

I'm working on an old archaeological field-survey dataset from the 1960's. I'm trying to create a polygon grid (fishnet) over a specific area where I can give the individual 'tiles'of the grid an address as attribute (for example grid square Aaa-1). I need this 'address' in order to link the specific squares with the data of the finds attached to this square (determinations, types, quantities, etc).

About the definitions of 'addresses': Finds are recorded in 1m2 squares within a grid of 100x100m meter. Along the x-axis the squares are numbered numerical 1-100 and along the y-axis alphanumerical (for example: Aaa-Aaj for the first 10 rows, Aba-Abj for the second 10 rows, etc until Ajj, which is nr 100.).

I've tried the FSC tools, with the OSGR Tool to create a grid, but this didn't worked out due to an inclination in my grid and it is not possible to position the grid on specific coordinates. The 1966 grid was not orientated to the true north, but set-out with an angle of 5,5 degrees in regards to the north. I have the coordinates of the four corners of the 100x100 meter square.

I'm currently working in QGIS, but if (easy) solutions are possible in other GIS environments, I'd be happy to hear.

  • Don't you just need to assign a grid-number, then have a number-name lookup table? – Vince Jul 7 '19 at 15:17
  • That's the approach i initially had in mind. I've tried the first step by a grid-cell-centroid-approach within pyqgis, however, due to the inclination of the grid, the grid numbers are not working (understandable of course) out as i hoped, i.e. the first row ends at grid-number 21 and 22 continues at row 2. On row 1 however, after 22, the numbering continues with 42. – Marten Jul 7 '19 at 18:04
  • A rotated grid is a pretty basic trigonometry problem. You might be trying for something too subtle. – Vince Jul 7 '19 at 18:18
0

I'm no programmer and there is probably a more elegant/efficient/'Pythonic' way of doing this than the code shown below. However, as a novice learning the ropes of PyQGIS, this was a fun task to undertake. I've tested this in QGIS version 2.18, if you need to convert it to work with the QGIS 3 API, hopefully someone on here can help you with that.

I define a function to transform the coordinates of each cell and create nested loops to populate a list of character strings.

Enter the coordinates of the south west corner in the 'coords' variable and all the points will be rotated around this point.

You add a new script to QGIS through toolbox > scripts > create new script, a full tutorial can be found here.

from qgis.core import *
from PyQt4.QtCore import QVariant
import string, math

#Parameters to set up
epsg = '27700' # Enter epsg code here - British National Grid shown as example ** Use a projected crs that uses metres for units
coords = [111111.11, 111111.11] # Enter the coordinates of the south west corner of the grid, this point will remain unchanged and all other points will be rotated around this
angleRadians = math.radians(-5.5) # Enter the rotation of the grid here in degrees. Enter a +ive value for CW rotation or -ive value for CCW rotation
maxCols = 100 # Maxiumum number of columns in the grid

#Create new layer and add a field called "ref"
lyrSetupStr = 'Polygon?crs=epsg:{}'.format(epsg)
lyrOut = QgsVectorLayer(lyrSetupStr, 'grid', 'memory')
lyrOut.startEditing()
newField = QgsField("ref", QVariant.String)
lyrOut.dataProvider().addAttributes([newField])
lyrOut.updateFields()

#function to transform points using trig
def func_rotate(x, y, degree):

    #this part works out transformation relative to start position 
    xOut = x * math.cos(angleRadians) - y * math.sin(degree)
    yOut = x * math.sin(angleRadians) + y * math.cos(degree)

    #return the transformation plus the coordinate of the starting grid to get the absolute position
    return xOut + coords[0], yOut + coords[1]

#create a list of lowercase ascii characters
alphaList = list(string.ascii_lowercase)

#create list to hold column references
rowIdList = []

#loop through first set of a-j rows
for i in range(10):

    #loop through second set of a-j rows and append to list
    for j in range(10):
        rowId = [alphaList[i], alphaList[j]]
        rowIdList.append(rowId)


#loop through rows
for rowNum in range(maxCols):
    #get the row reference from the list
    rowId = rowIdList[maxCols - rowNum - 1]

    #loop through columns
    for colNum in range(maxCols): 

        #call the transformation function for each of the 4 points of the cell square
        pt1 = func_rotate(maxCols - rowNum - 1, maxCols - colNum - 1, angleRadians)
        pt2 = func_rotate(maxCols - rowNum , maxCols - colNum - 1, angleRadians)
        pt3 = func_rotate(maxCols - rowNum , maxCols - colNum, angleRadians)
        pt4 = func_rotate(maxCols - rowNum -1, maxCols - colNum , angleRadians)

        #turn the returned coordinates into Well Known Text
        wkt = 'POLYGON (({} {}, {} {}, {} {}, {} {}))'.format(str(pt1[0]), str(pt1[1]), str(pt2[0]), str(pt2[1]), str(pt3[0]), str(pt3[1]), str(pt4[0]), str(pt4[1]))

        #Create a new feature, set the geometry 
        f = QgsFeature()
        geo = QgsGeometry.fromWkt(wkt)
        f.setGeometry(geo) 

        #add the ref field to the new feature and populate
        newFields = QgsFields()
        newFields.append(newField)
        f.setFields(newFields)
        attOut = '{}A{}{}'.format(str(maxCols - rowNum), str(rowIdList[colNum][0]), str(rowIdList[colNum][1]))
        f['ref'] = attOut

        #add the feature to the new layer
        lyrOut.dataProvider().addFeatures([f])

#add the new layer to the map
QgsMapLayerRegistry.instance().addMapLayer(lyrOut)
lyrOut.commitChanges() 
| improve this answer | |
  • Thank you for the script and the explanation, I'm eager to try the script. I'll let you know! – Marten Jul 11 '19 at 9:46

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.