I am trying to do a treatment on many feature files, so I would like to automate it.

In fact, I have one shapefile with the spatial distribution of some species, and one with vegetation type.

I would like to select (by attribute) one species in the Species Shapefile, then to select (by locality) all vegetation areas which intersect with their distribution area. Finally, I would like to have a shapefile with, as name, the species name, and with attributes and forms of the frequented vegetation types. And I would like to repeat this for all species (more than 100), and if possible do this in an easy manner (so it could be done by another person).

I have already tried this task using Sextante plugin, but I can't have the species name as a shapefile name at the end.

Can someone suggest a method for this?

  • 1
    From your description the entire work would be more suitable for a full-featured geodatabase like PostGIS or SpatiaLite. But a complete solution to perform what you need may be not trivial.
    – steko
    Commented Jan 26, 2013 at 19:58

2 Answers 2


This blog entry might help to understand how to do it in SEXTANTE:


Hope it helps

  • Hi, thanks for that, I want to try it, but I am not expert on this kind of script, so can you explain what we have to do with this script? where should we have to copy it? thanks.
    – Onesime
    Commented Jan 30, 2013 at 8:57
  • Great, thank you very much, I tried it (through the Python Console) and it works well. The next step, I will try to adapt it as a Sextante modeler. It is damage that there is no such command in such a tool as Sextante (set the name of the output file with some variables).
    – Onesime
    Commented Jan 31, 2013 at 13:46

This calls for a little script. To make it reproducible I would try to accomplish it in R. It should also be possible with QGis and Sextante by using the batch-execution (rightclick on function) in a Sextante model. Here you could first use the vector intersection tool and afterwards some kind of spatial join.

In R i would try it like this. You might have to modify the code as i don't know your data structure and variables.

library(raster);library(rgdal);library(sp)         # Load in required packages

vegetation <- readOGR("H:/Examplefolder",          # To load a vegetation polygon shape here 
                      "vegi")                      # called vegi.shp    

setwd(harddriveD)                                  # Now, switch to the directory containing your species shapes
                                                   # and use the following code 
species_files <- grep( ".shp",                     # to extract all shape names

                                                   # Now read in your speciesfiles in a loop 
for(name in species_files){                        # and do a  vegetation data
                                                   # overlay with your basename
    spec_name <- strsplit(name,split=".shp")[[1]]  # to get only the load in
                                                   # your species name shape. 

    spec_shp <- readOGR(".",spec_name)             # I assume that you have point data. Otherwise change the code.
    ov <- over(spec_shp,vegetation)                # Make a overlay with your vegetation data, 
                                                   # returns a dataframe of the overlaying vegetation shape attributes, 
                                                   # which are within your species shape. 
                                                   # This dataframe has the same number of rows as your input species shape. 
   cd <- coordinates(spec_shp);                    # Therefore you could just append the needed information to your species shape.
   proj <- proj4string(spec_shp)                   # save coordinates and proj.

                                                   # Append your vegetation data to your species shape
   spec_shp$Vegetation <- ov$ShrubSpecies          # Or whatever you want to get. 

   spp <- SpatialPointsDataFrame(                  # In the end export your shape again. 
                    coords=cd,                     # I would advise you to use a different folder. 
                    data=as.data.frame(spec_shp),  # In case you have polygons instead of Point data
                    proj4string=CRS(proj) )        # use the SpatialPolygonDataFrame function. -> help ?SpatialPolygonDataFrame
  writeOGR(spp,                                    #Should result in a new shape 
           "foldername",                           # which has your species name.
           driver="ESRI Shapefile")                      


I've made many assumption about your goal and the structure of your dataset. Most likely you have to correct the code to your needs before trying it out.

  • Thanks for your help, I will try it. My species data are in polygon (species distribution), but I think it's maybe quite the same? many thanks
    – Onesime
    Commented Jan 29, 2013 at 13:01
  • You just need to change some function (SpatialPolygonsDataFrame for instance) and most likely over returns a list of dataframes or something else.
    – Curlew
    Commented Jan 29, 2013 at 22:30

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.