I want to create a shapefile for each of the 72,000 census tracts in the United States that includes the census blocks in that tract using open source software.

I will probably start with a state shapefile that includes all of the blocks for that state and merge it with a GDB file. This will give me a list of blocks that are in the census tract.

But then I'm not sure how to split the state shapefile into smaller ones.

How can I do this?

Is there a way to automate this using open source software?

If so, what tools should I use?

I'm learning QGIS.

  • 1
    Does gis.stackexchange.com/questions/40798/… solve your problem? – BradHards Jun 17 '13 at 23:57
  • Maybe. I haven't used a programming language to work with shapefiles before. So I would need to learn how to read and write them - and what things like OGR are. I also don't know Python. I use PHP and MySQL. General method: load gdb and shapefile (in php if possible), merge them, parse it, find each block that is in the census tract, write out census tracts to a shapefile or kml. I'm planning on using kml as the final format. – Aaron Kreider Jun 18 '13 at 0:09
  • Can I open a geodatabase file in PHP? – Aaron Kreider Jun 18 '13 at 0:15
  • 2
    @AaronKreider PHP is for webpages, not scripting or programming. If you've learnt PHP, you can and should learn Python. There is a wide array of libraries available for python that you could use. No one builds GIS libraries for PHP because that's not what it's for - so you will have a hard time finding a GDB wrapper written in PHP. – nagytech Jun 18 '13 at 0:16
  • I found this library for processing shapefiles in PHP: github.com/fillerwriter/PHP5-Shape-File-Parser – Aaron Kreider Jun 18 '13 at 20:27

Basically, you can do the extract using ogr2ogr as long as you give the Census tract ID, so it's really an issue of getting 72,000 ogr2ogr calls.

ogr2ogr -where "tract = '<tract_id>'" /dest_folder /source_folder block_shapefile -nln block_shapefile_<tract_id>


  • You don't have to specify the source format, ogr2ogr will figure it out. You specified shapefile, so that's what I'm assuming.
  • If you are not changing the format, you also don't have to specify the destination format. Otherwise, to get a shapefile, add `-f "ESRI Shapefile"
  • I'm using the -where switch to subset the data. Remember that attribute query is much faster than spatial query. Your question was a little vague, so I don't know if you intend to join the blocks to the tracts by attribute or spatially, but I highly recommend the former.
  • Note that I am assuming your blocks shapefile has a tract ID column in it. Most Census data sets will have a hierarchy of ID columns, i.e. the blocks shapefile will probably have a state, county, and tract, as well as block ID column, but may not have a column concatenating them all together. They must be concatenated, because counties are only unique with states, tracts are only unique with counties, and blocks are only unique within tracts. So you will have to either (a) create a new field with state & county & tract as your ID field, or add all three criteria to the where clause.

So how do you build 72,000 ogr2ogr calls programmatically? You can use any tool you want, but here's an example with R:

dfBlocks = read.dbf("/source_folder/shapefile_name.dbf", as.is=TRUE)

strTract = unique(dfBlocks$tract_id)
for (i in 1:length(strTract)) {
  strOGR = paste(
    "ogr2ogr -where \"tract_id = '", strTract, 
    "'\" /dest_folder /source_folder layer_name -nln base_name_", 
    strTract, sep=""

You could also collect the system calls in one step, then iterate it to run the ogr2ogr call in a separate loop, or at a later time, or write it to a bash script that you run from the command line.

The major disadvantage is that it will scan the data source once for each ogr2ogr call, so actually importing the source data, iterating, and writing a shapefile for each row, would probably be more efficient. But I would recommend trying it in something other than R, which is somewhat slow at reading large spatial datasets (I ran an import just of the counties of the US, ~3000, and after 15 minutes I cancelled the import.)

| improve this answer | |
  • Hi Lee! I'am try use your script, but receive some troubles: I create in CentOS .sh file with content <code> #!/usr/bin/env Rscript library(foreign) dfBlocks = read.dbf("/source/farms.dbf", as.is=TRUE) strTract = unique(dfBlocks$name) for (i in length(strTract)) { strOGR = paste( "/usr/local/bin/ogr2ogr -overwrite -where \"name='",strTract,"'\" /destination/ /source/ farms -nln ",strTract,sep="" ) print (strOGR) system(strOGR) } <code>. When I run this script I receive only one first object from initial shapefile in destination folder. – spatialhast Oct 16 '14 at 11:43
  • Assume your shapefile has 1500 features, then strTract has 1500 feature names and length(strTract) is also 1500. for i in length(strTract) is saying for i in 1500, which is exactly one item (the last) rather than all items. Note in my code above the correct for statement is for i in 1:length(strTract) – Lee Hachadoorian Oct 18 '14 at 23:17
  • I avoid this troubles to add next code: <code> fileConn<-file("/rscript.sh") writeLines(c(strOGR), fileConn) close(fileConn) system('/rscript.sh')<code> – spatialhast Oct 20 '14 at 10: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.