2

I'm trying to build a grid over the world and cut it on an AOI. The initial grid have non trivial properties.

It is a grid of squares in EPSG:4326. The dimension is not round, they simply cut the world horizontally into 2048 squares ( => 1024 vertically) or .17578125° per cell. The point [0,0] being the lower left corner of a cell.

I wrote a Python script that manages to create this kind of grid :

from shapely.geometry import Point
import geopandas as gpd
import decimal as d
from itertools import product
import numpy as np

# the size is based on the planet grid size 
# the planet grid is composed of squared grid that split the world width in 2048 squares
diametre = 360/2048
radius = diametre/2
    
# compute the longitudes and latitudes for the whole world
longitudes = np.linspace(-180, 180, 2048)
latitudes = np.linspace(-90, 90, 1024)
    
# filter with the geometry bounds
min_lon, min_lat, max_lon, max_lat = aoi_gdf.total_bounds
longitudes = longitudes[(longitudes > (min_lon - radius)) & (longitudes < max_lon + radius)]
latitudes = latitudes[(latitudes > (min_lat - radius)) & (latitudes < max_lat + radius)]

# create the grid
points = []
for i, coords in enumerate(product(longitudes, latitudes)):
        
    x = d.Decimal(coords[0])
    y = d.Decimal(coords[1])
        
    points.append(Point(x, y))
    
# create a buffer grid in lat-long
grid = gpd.GeoDataFrame({'batch': batch, 'geometry':points}, crs='EPSG:4326') \
    .buffer(d.Decimal(diametre)) \
    .envelope \
    .intersection(aoi_shp_proj)
    
# filter empty geometries
grid = grid[np.invert(grid.is_empty)]

unfortunately when I watch my grid in QGIS I realized that the cells are actually overlaping. I also try without using the decimal lib but I end up with the same kind of problems (exact value of the Python float).

Is there a way to produce this kind of grid without creating ill shapes ?

EDIT

To show what I mean by ill shape. Here is the grid produce on Singapore : enter image description here

And here is a zoom on the corner

enter image description here

The tiles are in fact not aligned and overlap

1
  • can you show us a picture of what the problem is?
    – Ian Turton
    Commented Mar 3, 2021 at 17:01

1 Answer 1

0

Ok so I find a way to create the grid I wanted. 2 tricks where used:

  • First do not create buffer but construct the squares with the points coordinates to make sure that the same value is used (we cannot do arythmetic with floats unfortunately)
  • Second, in shp format the grid points where too precises and I was getting errors with the definitions of corners. GeoJson only use 13 decimals so the troncature was actually working. (reexport it to .shp afterward if needed)

NOTE: aoi_shp is a shapely geometry in EPSG:4326 projection. I added an extra step with a aoi_gdfto easily move from 1 projection to another

from pathlib import Path
from itertools import product

import geopandas as gpd
from shapely import geometry as sg
import numpy as np
from pyproj import CRS

# get the shape of the aoi in EPSG:4326 proj 
aoi_gdf = gpd.GeoDataFrame({'geometry': [aoi_shp]}, crs="EPSG:4326")
    
# retreive the bb 
aoi_bb = sg.box(*aoi_gdf.total_bounds)

# get the crs area of use 
crs_4326 = CRS.from_epsg(4326)
crs_min_x, crs_min_y, crs_max_x, crs_max_y = crs_4326.area_of_use.bounds

# create the set of point 
longitudes = np.linspace(crs_min_x, crs_max_x, 2048+1) # here you will have the points not the squares so +1
latitudes = np.linspace(crs_min_y, crs_max_y, 1024+1)

box_size = (crs_max_x-crs_min_x)/2048

# filter with the geometry bounds
min_lon, min_lat, max_lon, max_lat = aoi_gdf.total_bounds

# filter lon and lat 
lon_filter = longitudes[(longitudes > (min_lon - box_size)) & (longitudes < max_lon + box_size)]
lat_filter = latitudes[(latitudes > (min_lat - box_size)) & (latitudes < max_lat + box_size)]

squares = []
for coords in product(range(len(lon_filter)-1), range(len(lat_filter)-1)):
    
    ix = coords[0]
    iy = coords[1]

    squares.append(sg.box(
        lon_filter[ix], 
        lat_filter[iy], 
        lon_filter[ix+1], 
        lat_filter[iy+1]
    ))
        
# create a buffer grid in 4326 from the drawn squares
grid = gpd.GeoDataFrame({'geometry':squares}, crs='EPSG:4326')

# cut the grid to the country extends 
mask = grid.intersects(aoi_gdf.loc[0,'geometry'])
grid = grid.loc[mask]
    
# export the grid as a json file 
path = Path('~', 'tmp', 'test_grid.geojson').expanduser()
grid.to_file(path, driver='GeoJSON')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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