GDAL Grid – How to Create a Grid of 10degx10deg Polygons for a Shapefile Using GDAL

gdalpython

I want to create a shapefile using the GDAL library. And I want this Shapefile, when created, to have a grid of square polygons. And where each polygon cell is 10 degrees by 10 degrees, ultimately I want this grid to cover the entire world. As well, I also want the shapefile to contain a decimal column 'LAND', which should be populated with 0.0 for all rows/cells.

Having confusion here because I've never used GDAL before and just sifting through the documents is confusing me. I think I understand how to create a shapefile something like the following:

from osgeo import ogr
from osgeo import osr
import gdal

# specifying that I want to work with a shapefile
DriverName = 'My Shapefile'
driver = ogr.GetDriverByName(DriverName)


Filename = 'points.shp'

Just a little lost here. Not sure exactly how this would be done through GDAL.
I know if I am creating the world, the extent of the world is -180 deg to 180 deg longitude and -90 deg to 90 deg for latitude. And that I should create a nested loop somewhere to loop the outer of longitudes and the inner of latitudes, stepping a value of 10 each time for spatial resolution for the grid.

It's just the whole-set up and the library functionally is confusing not sure what way to go about all this.

The general structure of the code I'd think would be something like:
– Import GDAL and other modules
– create a spatial reference
– Create an empty shapefile
– create a layer in the shapefile and assign it the spatial reference
– create a set of point geometries as a ring
– add those point geometries in an order to a polygon geometry
– create a feature and add the polygon Geometry to it
– put the feature in the layer

I examine the GDAL library website and found this:
link: https://pcjericks.github.io/py-gdalogr-cookbook/

and then more specifically a section that said "Create a New Shapefile and Add Data"
link: https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#create-a-new-shapefile-and-add-data

I looked at the following but they did not help so much:
link: Creating square grid polygon shapefile with Python?

and

link: Creating grid polygons from coordinates using R or Python

I am trying to piece all this information together in order to solve my task but it is proving difficult.

Please help how would I create a shapefile being a grid of 10x10degree polygons?

Best Answer

This will do it (in R):

xmin <- -180
ymin <- -90
side <- 10 # degrees of the side of each cell
listP <- list() # creates an empty list (which is similar to an array)
lands <- vector() # creates an empty array
k <- 1
nx <- 36 # number of cols
ny <- 18 # number of rows
for (i in 1:nx) { # for each col
  x1 <- xmin + (i-1)*side
  x2 <- xmin + i*side
  print(paste(i,"/",nx,sep="")) # just in case it takes too long
  flush.console() # same above
  for (j in 1:ny) { # for each row
    y1 <- ymin + (j-1)*side
    y2 <- ymin + j*side
    p = Polygon(cbind(c(x1,x1,x2,x2,x1),c(y1,y2,y2,y1,y1)))
    pP = Polygons(list(p), as.character(k))
    listP[[length(listP)+1]] <- pP # append to the array
    lands[length(lands)+1] <- 0 # append to the array
    k <- k+1
  }
}
Grid <- SpatialPolygons(listP, 1:(k-1))
d <- data.frame(LAND=lands)
spdf <- SpatialPolygonsDataFrame(Grid, data=d)

Looks like the code in Java you showed should be something like:

form osgeo import gdal
from osgeo import ogr
from osgeo import osr
import os
import random

# Specify that you want to work with a Shapefile
DriverName = "ESRI Shapefile"
# Not sure what this line does
driver = ogr.GetDriverByName(DriverName)

# Name of shapefile to create
FileName = 'gride.shp'

# simple check to see if that file already exists in the system, if so delete it
if os.path.exists(FileName):
  driver.DeleteDataSource(FileName)

# Now create an Empty shapefile Data source
shapeData = driver.CreateDataSource(FileName)

fieldName = 'LAND'

lng_min = -180
lng_max = 180
step = 10

lat_min = -90
lng_max = 90

poly = org.Geometry(ogr.wkbPolygon)
for lng_var in range(lng_min, lng_max, step):
    for lat_var in range(lat_min, lat_max, step):
            ring = org.Geometry(ogr.wkbLinearRing)
            ring.AddPoint(lng_var, lat_var)
            ring.AddPoint(lng_var + step, lat_var)
            ring.AddPoint(lng_var + step, lat_var + step)
            ring.AddPoint(lng_var, lat_var + step)
            ring.AddPoint(lng_var, lat_var)
            poly.AddGeometry(ring)
Related Question