Main Page | Report Page

 

  Computers Forum Index » Computer Languages (IDL-PVWAVE) » Understanding IDLanROI...

Author Message
kBob...
Posted: Fri Oct 08, 2010 9:13 pm
 
I am attempting to use IDLanROI to determine if Shapefiles overlap
another Shapefile.

Sometimes it works, sometimes it doesn't. Especially, if the polygon
lies within the interior.

I put together a test program that retrieves the States and uses
IDLanROI to determine if it lies within the US.

Why does Colorado not considered to be in the United States, among
some other interior States?

Is IDLanROI not the IDL routine to use to perform this check?

Kelly Dean
Milliken, CO

;+
;
;
;
;----------------------------------------------------------------------
PRO ITT_PolyOverlap


states_path = FILEPATH('states.shp',SUBDIR=['resource', 'maps',
'shape' ])
country_path = FILEPATH('country.shp',SUBDIR=['resource', 'maps',
'shape' ])

PRINT, states_path
PRINT, country_path

oSHP = OBJ_NEW('IDLffShape', country_path )
ent1 = oSHP -> GetEntity( 34, /ATTRIBUTES ) ; USA
OBJ_DESTROY, oSHP

PRINT, ' -- Working with country : ', (*ent1.attributes).attribute_4

oSHP = OBJ_NEW('IDLffShape', states_path )
oSHP -> GetProperty, N_ENTITIES = num_ent
FOR staten = 0, num_ent DO BEGIN

oSHP = OBJ_NEW('IDLffShape', states_path )
ent0 = oSHP -> GetEntity( staten, /ATTRIBUTES )

oROI = OBJ_NEW('IDLanROI', (*ent1.vertices) )
pttest = oROI -> ContainsPoints( (*ent0.vertices) )
OBJ_DESTROY, oROI

overlap = WHERE( ptTest GT 0, count )

IF ( count GT 0 ) THEN BEGIN
PRINT, ' === ' , (*ent0.attributes).attribute_1, ' inside ',
(*ent1.attributes).attribute_4, ' === ', count
ENDIF ELSE BEGIN
PRINT, ' <<< ' , (*ent0.attributes).attribute_1, ' outside ',
(*ent1.attributes).attribute_4, ' >>> ', count
ENDELSE

ENDFOR

OBJ_DESTROY, oSHP

END
 
David Fanning...
Posted: Sat Oct 09, 2010 1:34 am
 
kBob writes:

Quote:
I am attempting to use IDLanROI to determine if Shapefiles overlap
another Shapefile.

Sometimes it works, sometimes it doesn't. Especially, if the polygon
lies within the interior.

I put together a test program that retrieves the States and uses
IDLanROI to determine if it lies within the US.

Why does Colorado not considered to be in the United States, among
some other interior States?

Is IDLanROI not the IDL routine to use to perform this check?

So, let's see, you are taking all the points from the Colorado
IDLanROI and asking the USA IDLanROI if any of those points are
exterior to the USA boundary, using the ContainsPoints method.
Is that right? And the answer is sometimes "yes"?

Just trying to get my head around it. :-)

Cheers,

David

--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
 
David Fanning...
Posted: Sat Oct 09, 2010 1:35 am
 
David Fanning writes:

Quote:
So, let's see, you are taking all the points from the Colorado
IDLanROI and asking the USA IDLanROI if any of those points are
exterior to the USA boundary, using the ContainsPoints method.
Is that right? And the answer is sometimes "yes"?

Oh, sorry. I saw your name at the "end" of the article,
and didn't realize there was code. I'll have another look.

David



--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
 
kBob...
Posted: Tue Oct 19, 2010 3:02 pm
 
On Oct 8, 3:13 pm, kBob <krd... at (no spam) gmail.com> wrote:
Quote:
 I am attempting to use IDLanROI to determine if Shapefiles overlap
another Shapefile.

 Sometimes it works, sometimes it doesn't. Especially, if the polygon
lies within the interior.

 I put together a test program that retrieves the States and uses
IDLanROI to determine if it lies within the US.

 Why does Colorado not considered to be in the United States, among
some other interior States?

 Is IDLanROI not the IDL routine to use to perform this check?

Kelly Dean
Milliken, CO

;+
;
;
;
;----------------------------------------------------------------------
PRO ITT_PolyOverlap

states_path = FILEPATH('states.shp',SUBDIR=['resource', 'maps',
'shape' ])
country_path = FILEPATH('country.shp',SUBDIR=['resource', 'maps',
'shape' ])

PRINT, states_path
PRINT, country_path

oSHP = OBJ_NEW('IDLffShape', country_path )
ent1 = oSHP -> GetEntity( 34, /ATTRIBUTES ) ; USA
OBJ_DESTROY, oSHP

PRINT, ' -- Working with country : ', (*ent1.attributes).attribute_4

oSHP = OBJ_NEW('IDLffShape', states_path )
oSHP -> GetProperty, N_ENTITIES = num_ent
FOR staten = 0, num_ent DO BEGIN

  oSHP = OBJ_NEW('IDLffShape', states_path )
  ent0 = oSHP -> GetEntity( staten, /ATTRIBUTES )

  oROI = OBJ_NEW('IDLanROI', (*ent1.vertices)   )
  pttest = oROI -> ContainsPoints( (*ent0.vertices) )
  OBJ_DESTROY, oROI

  overlap = WHERE( ptTest GT 0, count )

  IF ( count GT 0 ) THEN BEGIN
    PRINT, ' === ' , (*ent0.attributes).attribute_1, ' inside ',
(*ent1.attributes).attribute_4, ' === ', count
  ENDIF ELSE BEGIN
    PRINT, ' <<< ' , (*ent0.attributes).attribute_1, ' outside ',
(*ent1.attributes).attribute_4, ' >>> ', count
  ENDELSE

ENDFOR

OBJ_DESTROY, oSHP

END

To answer my own question ...

Yes, IDLanROI is NOT the IDL routine to perform this check.

After pondering on this for a couple of weeks, I figure what I needed
to do.

IDLgrROI is the prefer IDL routine to use.

It is a tessellation thing, especially when dealing with State and
Country Shapefiles that may contain gaps. You may get away with very
simple polygons using IDLanROI when dealing with a handful of
vertices, but IDLgrROI and IDLgrROIGroup are the routines you should
be using along with incorporating IDLgrTessellator in the code.

The code below produces the desire results I was after.

Kelly Dean
Milliken, CO

P.S. I thought IDL 8.0 removed memory leaks? I had to destroy some
Objects after a code run.

;+
;
; at (no spam) file_comments
; <P>Looking for (State) polygons that overlap (Country) polygons.
;
;----------------------------------------------------------------------
PRO ITT_TessOverlap

COMPILE_OPT DEFINT32, STRICTARR

states_path = FILEPATH('states.shp',SUBDIR=['resource', 'maps',
'shape' ])
country_path = FILEPATH('country.shp',SUBDIR=['resource', 'maps',
'shape' ])
oUSA = OBJ_NEW('IDLffShape', country_path )
oUSA -> GetProperty, N_ENTITIES = num_ent
;entUSA = oUSA -> GetEntity( 2, /ATTRIBUTES ) ; CAN
entUSA = oUSA -> GetEntity( 34, /ATTRIBUTES ) ; USA
;entUSA = oUSA -> GetEntity( 73, /ATTRIBUTES ) ; MEX
PRINT, ' -- Working with country : ', (*entUSA.attributes).attribute_4
oROIGroup = OBJ_NEW( 'IDLgrROIGroup' )
oROI = OBJARR( entUSA.n_parts )
FOR ii = 0L, entUSA.n_parts-1 do begin
IF (ii EQ entUSA.n_parts-1)THEN BEGIN
startpoint = (*entUSA.parts)[ii]
endpoint = entUSA.n_vertices-1
ENDIF ELSE BEGIN
startpoint = (*entUSA.parts)[ii]
endpoint = (*entUSA.parts)[ii+1]-1
ENDELSE
data = (*entUSA.vertices)[*, startpoint:endpoint]
oTess = OBJ_NEW('IDLgrTessellator')
oTess -> AddPolygon, data
x = oTess -> Tessellate(verts, polys)
oROI[ii] = OBJ_NEW('IDLgrROI', data = verts )
oROIGroup -> ADD, oROI[ii]
ENDFOR
oStates = OBJ_NEW('IDLffShape', states_path )
oStates -> GetProperty, N_ENTITIES = nstates
FOR staten = 0, nstates-1 DO BEGIN
overlap = 0
entState = oStates -> GetEntity( staten, /ATTRIBUTES )
PRINT, ' -- Working with state : ',
(*entState.attributes).attribute_1
FOR ii = 0L, entState.n_parts-1 do begin
IF (ii EQ entState.n_parts-1)THEN BEGIN
startpoint = (*entState.parts)[ii]
endpoint = entState.n_vertices-1
ENDIF ELSE BEGIN
startpoint = (*entState.parts)[ii]
endpoint = (*entState.parts)[ii+1]-1
ENDELSE
data = (*entState.vertices)[*, startpoint:endpoint]
oTess = OBJ_NEW('IDLgrTessellator')
oTess -> AddPolygon, data
x = oTess -> Tessellate(verts, polys)
pttest = oROIGroup -> ContainsPoints( verts )
npts = N_ELEMENTS( ptTest )
FOR sym = 1, 3 DO BEGIN
indx = WHERE( ptTest EQ sym, symcnt )
IF ( symcnt GT 0 ) THEN BEGIN
overlap = ( overlap GE 0 ) ? 1 : 0
ENDIF
ENDFOR
ENDFOR
IF ( overlap GT 0 ) THEN BEGIN
PRINT, ' === ' , (*entState.attributes).attribute_1, ' inside ',
(*entUSA.attributes).attribute_4, ' === '
ENDIF ELSE BEGIN
PRINT, ' <<< ' , (*entState.attributes).attribute_1, ' outside ',
(*entUSA.attributes).attribute_4, ' >>> '
ENDELSE
ENDFOR
OBJ_DESTROY, oROIGroup
FOR staten = 0, nstates-1 DO OBJ_DESTROY, oROI[staten]
PRINT, " -- That's all Folks!"
END
 
David Fanning...
Posted: Tue Oct 19, 2010 9:38 pm
 
kBob writes:

Quote:
To answer my own question ...

Yes, IDLanROI is NOT the IDL routine to perform this check.

After pondering on this for a couple of weeks, I figure what I needed
to do.

IDLgrROI is the prefer IDL routine to use.

It is a tessellation thing, especially when dealing with State and
Country Shapefiles that may contain gaps. You may get away with very
simple polygons using IDLanROI when dealing with a handful of
vertices, but IDLgrROI and IDLgrROIGroup are the routines you should
be using along with incorporating IDLgrTessellator in the code.

Can you elaborate a bit on why this "is a tessellation thing"?
I don't really understand what that means. And, supposing
I did understand what it means, what would it have to do
with choosing IDLgrROI over IDLanROI. Isn't IDLgrROI a
subclass of IDLanROI?

Cheers,

David


--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
 
Matt...
Posted: Mon Oct 25, 2010 4:17 pm
 
On Oct 22, 8:16 am, kBob <krd... at (no spam) gmail.com> wrote:

Quote:
To find the polygons that are holes, I use the IDL routine POLY_AREA.
The /SIGNED keyword returns the area as either positive or negative,
which tells you if the vertices are counterclockwise or clockwise. All
the State polygons from states.shp are clockwise.


Kelly,

Thanks so much for this tip. It has saved me. So other people out
there. Is this the preferred method to determine the vertex's
direction? internal vs external?

This is more out of interest than need. The poly_area is working like
a charm for me.

Thanks,
Matt

--
Matthew Savoie - Manager of Science Programmers
National Snow and Ice Data Center
(303) 735-0785 http://nsidc.org
 
 
Page 1 of 1    
All times are GMT
The time now is Thu Jul 31, 2014 5:38 pm