wtorek, 22 października 2013

Updating a field in a Geodatabase

I wanted to update in ArcGIS a point feature class with values representing polygons ID's in which the points are located:



I could do this by means of "Join Data" tool based on spatial location (ArcMap > Table of content > local menu > Joins and Relates > Join...). But this tool creates a new feature class which I do not need.

I could use "Spatial Join" tool from the ToolBox but then also a new feature class  is created ("The output must be unique from the input.").

Maybe there's a way in ArcGIS to just update like I wanted but I couldn't find it (some people say there isn't - see:

So I run SQL Server client and typed:

UPDATE Geobaza_08.dbo.PUNKT_ODWIERT
SET OBSZAR_ID = ob.Obszar_ID
FROM Geobaza_08.dbo.PUNKT_ODWIERT po, Geobaza_08.dbo.OBSZAR_BADAN ob
WHERE po.SHAPE.STIntersects(ob.SHAPE) = 1
;


I refreshed the table and I got what I needed:


I tried to make this in ArcGIS without using pure SQL. So I used Python and Arcpy (there should be an edit session opened before executing the following script):


#create update cursor 
uc = arcpy.da.UpdateCursor("PKTO",("OBSZAR_ID","SHAPE@","OBJECTID")) 
#iterate on points
for pkt in uc: 
     #select a polygon containing the point
     arcpy.SelectLayerByLocation_management("OBSZ","CONTAINS",pkt[1]) 
     #create select cursor from the selection
     obs = arcpy.SearchCursor("OBSZ")   
     #iterate on the select cursor 
     for ob in obs: 
          #write the selected polygon id into point field
          pkt[0] = ob.getValue("Obszar_ID")
          #update row in update cursor
          uc.updateRow(pkt)


Hmm... I rather prefer simple SQL update in the DBMS.

środa, 24 lipca 2013

Oracle - spatial join

Regarding spatial operations like SDO_ANYINTERACT and other you should remember: only one parameter can be a layer (spatially indexed geometry column), the other must represent one spatial object.

Let's try:


I want to know which objects from LAYER_1 has "anyinteract" realationship with the object 103 from LAYER_2.

select id
from layer_1
where SDO_ANYINTERACT(
        geom
        ,
        (select geom from layer_2 where id = 103)
        ) = 'TRUE'
;

        ID
----------
         1
         3

As soon as the second parameter represents one object we'll get the result. But if there are more objects (in this case we'll take two: object 102 and 103) – we won't. Let's try:

select id
from layer_1
where SDO_ANYINTERACT(
        geom
        ,
        (select geom from layer_2 where id in (103,102))
        ) = 'TRUE'
;


ORA-29902: error in executing ODCIIndexStart() routine
ORA-01427: single-row subquery returns more than one row
29902. 00000 -  "error in executing ODCIIndexStart() routine"
*Cause:    The execution of ODCIIndexStart routine caused an error.
*Action:   Examine the error messages produced by the indextype code and
           take appropriate action.

Proper way to do it is the spatial join as follows:

select x.id
from layer_1 x
  join layer_2 y on (SDO_ANYINTERACT(x.geom,y.geom) = 'TRUE')
where y.id in (103,102);


        ID
----------
         1
         3
         5

All the 'anyinteract' relationships we obtain by the following SQL:

select x.id LYR_1, y.id LYR_2
from layer_1 x
  join layer_2 y on (SDO_ANYINTERACT(x.geom,y.geom) = 'TRUE')
;

     LYR_1      LYR_2
---------- ----------
         1        103
         2        101
         3        103
         5        102

We can also display what exact type of relationship is this:

select
    x.id LYR_1
  , y.id LYR_2
  ,SDO_GEOM.RELATE(x.geom, 'determine', y.geom, 1e-2) as REL
from layer_1 x
  join layer_2 y on (SDO_ANYINTERACT(x.geom,y.geom) = 'TRUE')
;

    LYR_1      LYR_2         REL                      
---------- ----------------- -------------------------
         1        103        OVERLAPBDYINTERSECT 
         2        101        CONTAINS                
         3        103        OVERLAPBDYINTERSECT     
         5        102        OVERLAPBDYINTERSECT