[GRASS-dev] [GRASS GIS] #401: v.distance enhanced to report multiple categories

[GRASS GIS]

#401: v.distance enhanced to report multiple categories
-------------------------+--------------------------------------------------
 Reporter:  wilsonadam   |       Owner:  grass-dev at lists.osgeo.org
     Type:  enhancement  |      Status:  new                      
 Priority:  major        |   Milestone:  6.3.1                    
Component:  Vector       |     Version:  6.3.0 RCs                
 Keywords:  v.distance   |    Platform:  Linux                    
      Cpu:  x86-64       |  
-------------------------+--------------------------------------------------
 Summary:  I would like to request that the v.distance function be updated
 (with a flag?) to allow reporting of multiple categories in one layer.

 Example of why this is important:
 I am working with historical forest fire data and want to extract
 histories for different points.   I start with a shapefile that contains
 many (1000s) of polygons, many of which overlap in x-y space but vary over
 time (time information is in the attribute table). I could separate this
 into separate shapefiles (one for each year) prior to importing to GRASS,
 but this would result in almost 100 different layers and I would rather
 avoid it.  When I import it into grass using v.in.ogr, it is topologically
 cleaned and the result is a layer of (intersected) polygons, many of which
 have multiple categories that link to the attribute table.  For example, a
 single polygon could have burned in multiple years, so it is linked to
 multiple rows in the attribute table.  These multiple categories are
 visible with "v.category -g input=fire option=print" which results in
 something like this:
 2452

 2452/2540

 2452/2526/2540/2543

 2540/2575

 2406/2420/2517/2563/2581/2584

 Where each row is a unique polygon and the different elements are the
 various categories (rows in the attribute table) that are linked to it.
 So far so good.  But what I want to do is extract the fire history for a
 number of points, but v.distance only reports the last category for each
 polygon (which in my case is usually only the most recent fire) and
 reports "WARNING: more cats of to_layer."  So there seems to be a hidden
 ID value for each polygon (which would correspond to the invisible row
 number in the output above) but I cannot seem to access it directly.  If I
 could, then it would be possible to v.distance to that ID, then use the
 output above to link a given point to several fire records.

 If v.distance was updated to include multiple categories in the same
 layer, I would be able to do this easily. This has been proposed before:
 http://www.mail-archive.com/grass-user@lists.osgeo.org/msg02056.html.  I
 would like to encourage this revision (though maybe with a -m flag so you
 could turn this feature on if wanted).  It would ideally (for me) return a
 table with multiple records for each point, each with a different category
 from the polygon layer.  For example, something like this:

 point | fireyear

 1           1950

 1           1975

 1           2002

 2           1960

 2           1972

 3           1954

 For reference, I have been able to do this quite easily with a loop in R
 with the following code:
 (though this is specific to my dataset, some changes will be needed).

 ###################  This code intersects a set of points with any number
 of polygons (which may be overlapping)

 library(sp);library(rgdal)
 fires=readOGR("/media/Data/Work/Regional/CFR/FireAnalysis/FireData/fires_15072008","all_fires_07_08")
 #read in shapefile with overlapping polygons
 d=as.data.frame(cbind(slot(fires,"data")[1,],point=1))  #use first row as
 template, add "point" as placeholder to be filled later
 nfires=nrow(slot(fires,"data")) #get the number of polygons
 for(i in 1:nfires) {  #loop through each polygon one at a time
   d2=overlay(fires[i,],points)  #do overlay of all points for each fire
 polygon (this may result in lots of NAs)
   d2$point=as.factor((1:nrow(d2))+100)  #add point ID - my point IDs start
 at 101 and go up, you will have to adjust this
   d=rbind(d,d2)  #bind this polygon's overlay to the previous one
   print(paste(i," out of ",nfires))  #print progress
 }

 d=d[-1,] #remove first line - used to start dataframe
 d=d[!is.na(d$FIREREFERE),] #get rid of all the NAs using a field that is
 always populated
 d2=merge(d,slot(points,"data"),by.x="point",by.y="Locality_n",all.x=T)
 #merge with point data to get point attributes for each point

 #########################################

-- 
Ticket URL: <http://trac.osgeo.org/grass/ticket/401>
GRASS GIS <http://grass.osgeo.org>