0

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?

1

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)
  • Thanks! Could you explain a bit what is happening here? How would a 10deg x 10deg polygons be forged from this script? Should this be at some point outputting a final shapefile. I should to be able to open it up on QGIS. – Ki_ Apr 26 '18 at 16:19
  • What are the <- doing exactly, sorry if that's a dumb question. Could I use '=' instead? – Ki_ Apr 26 '18 at 16:20
  • @Ki_ They are almost the same, but = is more used in narrower scopes, like in function parameters. I don't know the exact difference, thus I always use <-, except in function definitions. And no, it's not a dumb question at all. – Rodrigo Apr 26 '18 at 16:23
  • @Ki_ The loop is creating an array of Polygons (listaP) composed of one or more Polygon (p) (one, in this case), that will later be assembled together in the last 3 functions. If you look at the slotNames of any shapfile, you'll get closer to understand how they work. – Rodrigo Apr 26 '18 at 16:26
  • I am not sure what tweaks to make here I will show what it should output, these were the directions: imgur.com/a/nv3n8Kr – Ki_ Apr 26 '18 at 16:41

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.