Skip to content

converting bounding box to shapefile problem

7 messages · Barry Rowlingson, Roger Bivand, mihai.niculita +1 more

Original
mihai.niculita
# 
Hy all,

I need to convert the bounding box of several rasters to shapefile. I wrote
a script starting from several scripts available over the Internet, and I am
stuck at the final command which should write the bounding box polygons
(contained in spatpoldf.list) obtained with sp package to shapefile using
rgdal:
+ write.list[[i]] <- writeOGR(spatpoldf.list[[i]], 
+ ".", layer=map.list[[i]], driver="ESRI Shapefile")
+ }
Error in writeOGR(spatpoldf.list[[i]], ".", layer = map.list[[1]], driver =
"ESRI Shapefile") : 
  number of objects mismatch
map.list is a vector containing the raster names and looks like this:

[1] "N21W157.asc" "N21W158.asc" .....

I mention that the bounding box polygon does not have CRS because I can not
attach decimal degree coordinates with sp package. I really need the
Lat/Long coordinates because finally I must apply the script to a big number
of SRTM1 rasters which cover several UTM zones.

The entire script and the data is available at:
http://www.geomorphologyonline.com/script.r
http://www.geomorphologyonline.com/data.tgz

I work with R.2.12.1 x64 on Windows and I have the latest sp and rgdal for
x64.

Thank you in advance for any help with this!

--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/converting-bounding-box-to-shapefile-problem-tp6462495p6462495.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
Barry Rowlingson
#
On Fri, Jun 10, 2011 at 3:38 PM, mihai.niculita <mihai.niculita at uaic.ro> wrote:
This is because the SpatialPolygonDataFrames you try to write to
shapefiles have 7 items in their data but only one polygon - a square.
The data should only be a data frame with one row.

 It comes from this:

ID <- array(map.list, dim=c(7,1))   # why 7 here?
IDD <- as.data.frame(ID)

and then:

spatpoldf.list[[i]] <- SpatialPolygonsDataFrame(spatpol.list[[i]],
IDD, match.ID = FALSE)

I'm not sure why SpatialPolygonsDataFrame doesnt complain at this
point. I do know that if you replace IDD with data.frame(ID=1) or some
other single-row data frame when you construct it then it will work.

[Also, you only need about three of those packages loaded in your file]

Barry
Roger Bivand
#
On Fri, 10 Jun 2011, Barry Rowlingson wrote:

            

      
      
      
    

        
It doesn't complain in principle because the user turned off ID checking, 
but you are right that nrow(df) should be checked even if ID checking is 
turned off.

Thanks for noticing an open barn door!

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      mihai.niculita
    

    

      
      #
    

  


  

      
Thank you Barry and Roger,

You are right with the objects. This part from sp package I didn't
understood to well. My code is trying to obtain a list of
SpatialPolygonsDataFrame:

ID <- array(map.list, dim=c(7,1))
IDD <- as.data.frame(ID)
spatpol.list <- list()
polygon.list <- list()
spatpoldf.list <- list()
for(i in 1:length(bbx.list)) {
polygon.list[[i]] <- list(Polygons(list(Polygon(bbx.list[[i]])),
map.list[[i]]))
spatpol.list[[i]] <- SpatialPolygons(polygon.list[[i]])
spatpoldf.list[[i]] <- SpatialPolygonsDataFrame(spatpol.list[[i]], IDD,
match.ID = FALSE)
}

I mention that I want to process automatically almost 800 rasters. In this
case I had 7 rastesr, that's why ID is an array of 7 rows.
The sp documentation says for the above code:
"srl       list with Polygon-class objects" == which in my script is
bbx.list[[i]]
"ID       character vector of length one with identifier" == !!!!this I
didn't understood, and in my case I choose the vector with the names of arc
ascii grid files
"Srl       list with objects of class Polygons-class" == which in my script
is polygon.list[[i]]
"pO       integer vector; plotting order; if missing in reverse order of
Polygons area" 
"Sr        object of class SpatialPolygons-class" == which in my script is
spatpol.list[[i]]
"data     object of class data.frame; the number of rows in data should
equal the"
            number of Polygons-class objects in Sr" ==!!!!this I didn't
understood, and in my case I transformed the vector with the names of arc
ascii grid files in a matrix
"match.ID       logical: (default TRUE): match SpatialPolygons member
Polygons ID slot values"
                    with data frame row names, and re-order the data frame
rows if necessary."
                    If character: indicates the column in data with Polygons
IDs to match"

In my case spatpol.list[[i]] should have i spatial polygons, not only one as
Barry say.
As a conclusion if some one can give me an example of how "ID" and "data"
should look I think I will resolve the problem.


-----
teaching assistant/phD student
Department of Geography
Faculty of Geography and Geology
University Al. I. Cuza
Iasi, Romania
email: mihai.niculita at uaic.ro
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/converting-bounding-box-to-shapefile-problem-tp6462495p6463373.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

  
  


  


      

    

    

      
      

        


  

        

      Barry Rowlingson
    

    

      
      #
    

  


  

      
On Fri, Jun 10, 2011 at 7:34 PM, mihai.niculita <mihai.niculita at uaic.ro> wrote:

        

            

      
      
      
    

        
Well, no, spatpol.list is a list of the same length as bbx.list. Each
element of spatpol.list, referred to as spatpol.list[[i]] is a single
spatial polygon (with four coordinates).

            

      
      
      
    

        
You have to decide if you want to create 7 different shapefiles, one
for each of your grids and each containing one bounding box polygon,
or one shapefile with 7 features.

 You are very close to the first solution. A shapefile, and a
SpatialPolygonsDataFrame, contains 1 or more features, where each
feature has a geometry (the polygon) and some attributes (the data
frame). Currently you create your SpatialPOlygonsDataFrame from a
single geometry and a data frame (IDD) with 7 rows - these dont match.
You need to make sure that IDD is only *one* row.


Barry

  
  


  


      

    

    

      
      

        


  

        

      Robert J. Hijmans
    

    

      
      #
    

  


  

      
Interesting problem

            

      
      
      
    

        
Here is a solution, I think, if you want a single shapefile with all
bounding boxes

library(raster)
maps <- list.files(, ".asc", full.names=TRUE)

# get extents. Advantage of this no need to read the raster file values
ext <- lapply(maps, function(x)  extent(raster(x)))

# get SpatialPolygons (is there an easier way to append SpatialPolygons with
the same ID?)
sp <- SpatialPolygons(lapply(1:length(ext), function(x) { Polygons(list(
as(ext[[x]], 'SpatialPolygons')@polygons[[1]]@Polygons[[1]]), x) } ))

# data.frame
dat <- do.call(rbind, lapply(ext, function(x)  as.vector(bbox(x))))
dat <- data.frame(dat, filename=basename(maps))
colnames(dat)[1:4] = c('xmin', 'ymin', 'xmax', 'ymax')

sp <- SpatialPolygonsDataFrame(sp, dat)

writeOGR(sp, ".", layer='filename', driver="ESRI Shapefile")

plot(sp)
xy = coordinates(sp)
text(xy[,1], xy[,2], 1:nrow(xy))


--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/converting-bounding-box-to-shapefile-problem-tp6462495p6464586.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

  
  


  


      

    

    

      
      

        


  

        

      mihai.niculita
    

    

      
      #
    

  


  

      
Thanks Barry and Robert,

Now I understood the meaning of the ID and IDD refering to
SpatialPolygonsDataFrame function. And I managed to run the entire script.
My intention is to obtain a shapefile of the bounding box for every raster
file. I think I can modify Robert shortcut to obtain this also.

Thanks again for your time!

-----
teaching assistant/phD student
Department of Geography
Faculty of Geography and Geology
University Al. I. Cuza
Iasi, Romania
email: mihai.niculita at uaic.ro
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/converting-bounding-box-to-shapefile-problem-tp6462495p6464692.html
Sent from the R-sig-geo mailing list archive at Nabble.com.