3

I have Baseline points. Around these points I would like to generate a grid that is a at a specified distance. Lets say that this grid will be 10km x 10km. It should have a limited amount of rows and columns for-instance 3 x 3. Please see what I mean in the attached picture.

I can generate a Fishnet and then move the grids by hand but that will be labour intensive due to the amount of Baseline points I have.

enter image description here

2
  • In the pic the baseline point is not in the center of the 3x3-grid. Which is the requirement concerning baseline point position and grid? Do other points play a role? How many grids are to be created per point layer? Are they allowed to overlap, in case you want many points surrounded by grid cells? Should cells be attributed in order to relate them to "their" baseline point? Is creation of surrounding center points optional/mandatory, or an alternative to the grid?
    – Detlev
    Mar 5, 2016 at 16:29
  • It should be in the centre of the 3x3 grid. Its just a quick illustration. Other points don't play a role. Each baseline point is on its own and even with the grid the supporting points wont overlap each other. And if they do its not a problem. The grid and amount of points is exactly as displayed. 3x3 grid with the baseline point forming the centre and 8 points surrounding it. Yes it will be nice if the supporting points can take up a unique attribute from the baseline point. The end result should be Points so yes I would say that the points is mandatory, Grid just to display the idea Mar 6, 2016 at 5:08

2 Answers 2

1

In QGIS, you can set up a custom CRS with a transversal mercator projection based on lon_0 and lat_0 of your base point:

+proj=tmerc +lon_0=-168 +lat_0=59 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

Change the project CRS to that custom CRS, and creat a vector grid from -15000 to +15000 in both directions. Make sure the grid is added to the legend with your custom CRS. Then densify the grid by 29, and reproject it to your UTM zone or whatever project CRS you are usually working with.

For other base points, you just have to edit the lon_0 and lat_0 of your custom CRS. You can still use the same grid file. You only need the WGS84 degree coordinates of your base points.

1

Based on Andre's suggestion to transform each base point to TM projection with itself as the origin, you can automate the generation of the 3x3 grid points with the following Python script.

Each grid point gets the feature id of its base point. Matrix definition and distance to base points are defined as variables, so feel free to change them.

Copy&paste into Python console or python editor window, and run.

# get base point layer
bases = QgsMapLayerRegistry.instance().mapLayersByName('Point')[0]

# define result layer
result = QgsVectorLayer('Point?crs=EPSG:4326&field=id:integer', '3x3Grid', 'memory')
prov = result.dataProvider()
fields = prov.fields()

# define the 3x3 matrix; others are possible as well
dist = 10000  # meters
mat = [[(-dist, dist), (0, dist), (dist, dist)], 
       [(-dist, 0), (0, 0), (dist, 0)], 
       [(-dist, -dist), (0, -dist), (dist, -dist)]]

# set target coordinate reference system
trs = QgsCoordinateReferenceSystem()

feats = []
for base in bases.getFeatures():
    # get coordinates of projection origin
    lon_0 = base.geometry().asPoint().x()
    lat_0 = base.geometry().asPoint().y()
    # get copy of basepoint
    geom = QgsGeometry(base.geometry())
    # set projection parameters
    trs.createFromProj4('+proj=tmerc +lon_0=%f +lat_0=%f +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs' % (lon_0, lat_0)) 
    transformation = QgsCoordinateTransform(bases.crs(), trs)

    # transform point to target coordinate system
    geom = transformation.transform(lon_0, lat_0, direction=QgsCoordinateTransform.ForwardTransform)

    # create coordinates of matrix points
    x_0 = geom.x()
    y_0 = geom.y()
    coords = [(x_0+mat[r][c][0], y_0+mat[r][c][1]) 
              for r in range(len(mat)) 
              for c in range(len(mat[r]))]

    for coord in coords:
        feat = QgsFeature(fields)
        # transform matrix points back to original projection
        point = transformation.transform(coord[0], coord[1], direction=QgsCoordinateTransform.ReverseTransform )
        feat.setGeometry(QgsGeometry().fromPoint(point))

        # set id to id of base point
        feat['id'] = base.id()
        feats.append(feat)

# add all feats
prov.addFeatures(feats)
# add result to legend
QgsMapLayerRegistry.instance().addMapLayer(result)

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.