1

I'm working on a script where I create a grid with the bounding box of a geodataframe, a width and a height (they might be equals for a square grid).

I don't know why but when I try to make a 10m square grid my script generate a 7m square grid.

All my process works with geodataframe in Pseudo Mercator (3857) who have meter for unit.

Here is the piece of script for the grid:

Xmin, Ymin, Xmax, Ymax = buffer_CH.total_bounds # coordonées de la bounding box du buffer de l'enveloppe concave

    cols = abs(int((Xmax - Xmin) / W))+1 # nombre de colonnes
    rows = abs(int((Ymin - Ymax) / H))+1 # nombre de lignes

    XleftOri, XrightOri, Ybottom, Ytop = Xmin, Xmin+W, Ymin, Ymin+H
    polygons = []
    print(XleftOri, XrightOri, Ybottom, Ytop)
    for r in range(rows):
        Xleft = XleftOri
        Xright = XrightOri
        for c in range(cols):
            polygon = Polygon([(Xleft,Ybottom),(Xright,Ybottom),(Xright,Ytop),(Xleft,Ytop)])
            polygons.append(polygon)
            Xleft = Xright
            Xright += W
        Ybottom = Ytop
        Ytop += H

    grid = gpd.GeoDataFrame({'geometry':polygons})
    grid.crs = {'init': 'epsg:3857'}
    grid.to_file("{0}/{1}_{2}_zonage-carre/{3}-grid.shp".format(os.chdir, date, folder_name, shp_name))

Issue seems to come from this operation:

rotate_grid = grid.rotate(angleinverse, origin=(0,0), use_radians=False) # Rotation
rotate_grid = gpd.GeoDataFrame(rotate_grid) # Conversion de la GeoSerie en GeoDataFrame
rotate_grid = rotate_grid.rename(columns={0:'geometry'}).set_geometry('geometry') # Identification de la colonne geometry
rotate_grid.crs = {'init': 'epsg:3857'} # attribution du crs
rotate_grid.to_file("{0}/{1}_{2}_zonage-carre/{3}-RotateGrid.shp".format(os.chdir, date, folder_name, shp_name))
  • 1
    how are you measuring the grid size? – Ian Turton Jul 30 '18 at 15:34
  • The measurement tool is likely returning the 'geodesic' distance which is about cos(latitude) shorter east-west than the Mercator distance. Do you get close to 10 m if you measure north-south? – mkennedy Jul 30 '18 at 17:12
  • I measure the grid size with QGIS measurement tool. I'm closer to 10m when I measure east-west but just 1 or 2cm larger. – Tim C. Jul 31 '18 at 7:27
1

Before saving to file assign the desired crs to your GeoDataFrame:

from fiona.crs import from_epsg
grid = gpd.GeoDataFrame({'geometry':polygons})
grid.crs = from_epsg(3857)

EDIT:

instead of using qgis measurement tool to compute the lenth of your polygons, use compute the length of the exterior Linear ring:

rotated_grid['perimetre'] = rotated_grid['geometry'].apply(lambda x: x.exterior.length)

so you can check if indeed it the polygons have changed size, or is it just an error due to the qgis measurement tool

  • Yes I forgot that step. Effectively it's better, closer to 10m, 9.82m. But grid export isn't the final of my process and I do many things on it after: rotation, spatial join, 1st clip, 2nd clip. And at the end, my clipped grid is already wrong sized... – Tim C. Jul 31 '18 at 7:42
  • you can apply that CRS to buffer_CH and every GeoDataFrame on your code, to make sure your coordinates won't go wrong – Hicham Zouarhi Jul 31 '18 at 7:51
  • I tried it, it didn't change anything. So I verified my process step by step and squares size changed during the grid rotation. I'll add this code part in my main post. – Tim C. Jul 31 '18 at 8:12
  • @TimC. are you sure that angleinverse is in degrees ? – Hicham Zouarhi Jul 31 '18 at 8:27
  • Yes, I verified it, it's ok. – Tim C. Jul 31 '18 at 8:34

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.