Dear all,
currently I am developing an pure IDL (No ENVI functions) based remote
sensing application which involves supervised classification.
Therefore I need to import a Geotiff with the raster data plus an ESRI
shapefile (Polygons) with the trainings-data. I then converted the
shape file to an IDLanROI object and build a mask using
IDLanROI::ComputeMask method. At this point I got stuck as I don't
know how spatially link both image and roi even when they are in UTM
projection.

Here is my code so far:
;Import Multispectral Geo TIFF
img=READ_TIFF('some_multispectral.tiff',GEOTIFF=geokeys,interleave=2,)
s= SIZE (img, /DIMENSIONS)

;Setup Map projection based on GeoTIFF tags  (UTM 17N, Wgs84)
mapCoord = GeoCoord(img(*,*,1),geokeys)
mapStruct = mapCoord -> GetMapStructure()

;import shapfile (aggain UTM 17N, Wgs84)
myshape=OBJ_NEW('IDLffShape', 'some_esri_multipolygone.shp')

;Get the first polygon
ent=myshape->IDLffShape::GetEntity(1)

;Convert SHP to ROI object
myroi = OBJ_NEW( 'IDLanROI',( *ent.vertices)[0,*] ,(*ent.vertices)
[1,*] )

;Apply ComputMask method
maskResult = myroi -> ComputeMask(DIMENSIONS = [s[0],
s[1]],MASK_RULE=2)

; Here starts the problem;  I have a 5000x5000 raster (img) and the
image coordinate system start with [0,0] but the myroi coordinates are
much higher (eg  227984.00 ; 991472.00) Thus everything is masked.
;There is a keyword LOCATION for ComputeMask but I am not sure how to
use it and where to get the right values for it

index= WHERE(maskResult gt 0)
subset=img(index)

OBJ_DESTROY, myshape
OBJ_DESTROY, myroi

Paul Magdon writes: 

> ;Convert SHP to ROI object
> myroi = OBJ_NEW( 'IDLanROI',( *ent.vertices)[0,*] ,(*ent.vertices)
> [1,*] )
> 
> ;Apply ComputMask method
> maskResult = myroi -> ComputeMask(DIMENSIONS = [s[0],
> s[1]],MASK_RULE=2)
> 
> ; Here starts the problem;  I have a 5000x5000 raster (img) and the
> image coordinate system start with [0,0] but the myroi coordinates are
> much higher (eg  227984.00 ; 991472.00) Thus everything is masked.
> ;There is a keyword LOCATION for ComputeMask but I am not sure how to
> use it and where to get the right values for it

I have to admit I have never done this before. (Although
I would like to try if you could make the image and 
shape file available to me.) But, it seems to me
that when you create the ROI, you will have to give
the ROI the range of the image to make this work.

   mapCoord -> GetProperty, XRANGE=xr, YRANGE=yr
   myroi = OBJ_NEW( 'IDLanROI',( *ent.vertices)[0,*] ,$
      (*ent.vertices)[1,*], ROI_XRANGE=xr, $
      ROI_YRANGE=yr )

That way, the ROI and the image are mapped to the same
coordinate system.

I'm curious to know if that works. :-)

Cheers,

David

P.S. The next book I am thinking of writing is a small
book about map projections in IDL. The more examples I
have like this, the better. I have to write this book
before I forget everything I've learned in the past four
years. :-)


-- 
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.")

On 23 Dez., 22:34, David Fanning <n...@dfanning.com> wrote:
> Paul Magdon writes:
> > ;Convert SHP to ROI object
> > myroi =3D OBJ_NEW( 'IDLanROI',( *ent.vertices)[0,*] ,(*ent.vertices)
> > [1,*] )
>
> > ;Apply ComputMask method
> > maskResult =3D myroi -> ComputeMask(DIMENSIONS =3D [s[0],
> > s[1]],MASK_RULE=3D2)
>
> > ; Here starts the problem; =A0I have a 5000x5000 raster (img) and the
> > image coordinate system start with [0,0] but the myroi coordinates are
> > much higher (eg =A0227984.00 ; 991472.00) Thus everything is masked.
> > ;There is a keyword LOCATION for ComputeMask but I am not sure how to
> > use it and where to get the right values for it
>
> I have to admit I have never done this before. (Although
> I would like to try if you could make the image and
> shape file available to me.) But, it seems to me
> that when you create the ROI, you will have to give
> the ROI the range of the image to make this work.
>
> =A0 =A0mapCoord -> GetProperty, XRANGE=3Dxr, YRANGE=3Dyr
> =A0 =A0myroi =3D OBJ_NEW( 'IDLanROI',( *ent.vertices)[0,*] ,$
> =A0 =A0 =A0 (*ent.vertices)[1,*], ROI_XRANGE=3Dxr, $
> =A0 =A0 =A0 ROI_YRANGE=3Dyr )
>
> That way, the ROI and the image are mapped to the same
> coordinate system.
>
> I'm curious to know if that works. :-)
>
> Cheers,
>
> David
>
> P.S. The next book I am thinking of writing is a small
> book about map projections in IDL. The more examples I
> have like this, the better. I have to write this book
> before I forget everything I've learned in the past four
> years. :-)
>
> --
> 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.")


Dear David,
sorry for  my late reply but as we have Christmas there are other
duties as well :) . By the way: merry Christmas!
I tried your suggestion to use ROI_XRANGE and got the following Error:
% Keyword ROI_XRANGE not allowed in call to: OBJ_NEW
As fare as I understand the IDL documentation  ROI_XRANGE can not be
set for an IDLanROI object.
Today I will try to use the LOCATION keyword or find an other option
to transfer the UTM coordinates into proper pixel coordinates.
I will send you a shp + corresponding raster as soon as possible.
Cheers,
Paul

Paul Magdon writes: 

> I tried your suggestion to use ROI_XRANGE and got the following Error:
> % Keyword ROI_XRANGE not allowed in call to: OBJ_NEW
> As fare as I understand the IDL documentation  ROI_XRANGE can not be
> set for an IDLanROI object.

Ah, right, I didn't notice that. Sorry. :-)

OK, then I suppose you have to convert the pixel
locations into the same units as the ROI (probably
lat/lon). If your image coordinates are in projected
XY space, then you put them into lat/lon space with
Map_Proj_Inverse. Map_Proj_Forward goes the other way.

I tend to *always* work in projection XY space, so
I would probably convert the shape file to these
coordinates before I stored them in the ROI.

Anyway, if you send me the files, I'll write up an
example for you, although I can't promise to be
quick about it, unfortunately. Too much going on 
here, and money is running out faster than my
fingers can type. :-)

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.")

On Dec 25 2010, 4:34=A0pm, David Fanning <n...@dfanning.com> wrote:
> Paul Magdon writes:
> > I tried your suggestion to use ROI_XRANGE and got the following Error:
> > % Keyword ROI_XRANGE not allowed in call to: OBJ_NEW
> > As fare as I understand the IDL documentation =A0ROI_XRANGE can not be
> > set for an IDLanROI object.
>
> Ah, right, I didn't notice that. Sorry. :-)
>
> OK, then I suppose you have to convert the pixel
> locations into the same units as the ROI (probably
> lat/lon). If your image coordinates are in projected
> XY space, then you put them into lat/lon space with
> Map_Proj_Inverse. Map_Proj_Forward goes the other way.
>
> I tend to *always* work in projection XY space, so
> I would probably convert the shape file to these
> coordinates before I stored them in the ROI.
>
> Anyway, if you send me the files, I'll write up an
> example for you, although I can't promise to be
> quick about it, unfortunately. Too much going on
> here, and money is running out faster than my
> fingers can type. :-)
>
> 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.")

Hi,
In my case both the shape and the raster are in UTM projection which
is already a Cartesian projection. So why should I convert them into a
lat/long projection?
@David: I send send you the example files or is there an option to
'post' them here?
Cheers Paul

Paul Magdon writes: 

> In my case both the shape and the raster are in UTM projection which
> is already a Cartesian projection. So why should I convert them into a
> lat/long projection?

Right. Stay away from lat/lon projections if you want
things to make sense. :-)

> @David: I send send you the example files or is there an option to
> 'post' them here?

You can just send them to me or put them in a place
where I can pick them up. I am reluctant to publish
my e-mail on a newsgroup, since my spam will increase
10-fold. But my name is David and I live at the same
place as Coyote. :-)

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.")

On 1 Jan., 19:32, David Fanning <n...@dfanning.com> wrote:
> Paul Magdon writes:
> > In my case both the shape and the raster are in UTM projection which
> > is already a Cartesian projection. So why should I convert them into a
> > lat/long projection?
>
> Right. Stay away from lat/lon projections if you want
> things to make sense. :-)
>
> > @David: I send send you the example files or is there an option to
> > 'post' them here?
>
> You can just send them to me or put them in a place
> where I can pick them up. I am reluctant to publish
> my e-mail on a newsgroup, since my spam will increase
> 10-fold. But my name is David and I live at the same
> place as Coyote. :-)
>
> 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.")

Hi David,
did you received my mail with the test dataset?  I don't want to
stress you it's just that I am not sure if I used the right E-mail
adress :)
Cheers
 Paul

Paul Magdon writes: 

> did you received my mail with the test dataset?  I don't want to
> stress you it's just that I am not sure if I used the right E-mail
> adress :)

Yes, sorry, I received it. Just pretty busy right now. 
I'll get to 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.")

Hi Paul,

If what you want is a mask where the pixels within each polygon have
as DN its ID, there is a simpler way which doesn't require the
IDLanROI class:

;Import GeoTIFF and create mask array
img=READ_TIFF('some_multispectral.tif',GEOTIFF=geokeys,interleave=2)
psz= mapinfo.ModelPixelSCALETAG[0] ; pixel size
x0= geokeys.ModelTiePointTag[3] - geokeys.ModelTiePointTag[0]*psz
y0= geokeys.ModelTiePointTag[4] + geokeys.ModelTiePointTag[1]*psz
; NB, x0 and y0 are respectively the easting and northing
; of the NW corner of the image
s= SIZE(img, /DIMENSIONS)
mask= LONARR(s)

;import shapefile
myshape= OBJ_NEW('IDLffShape', 'some_esri_multipolygone.shp')

; populate the mask (assumes is in the same projection
;and covers the same extent as the geotiff)
myshape->IDLffShape::GetProperty, N_ENTITIES=n
FOR i=0L, n-1 DO BEGIN
  feati=  myshape->IDLffShape::GetEntity(i)
  featix= Round((Reform((*feati.vertices)[0,*])-x0)/psz)
  featiy= Round((y0-Reform((*feati.vertices)[1,*]))/psz)
  featis= POLYFILLV(featix, featiy, ns, nl)
  IF featis[0] NE -1 THEN mask[featis]= feati.ishape +1
ENDFOR

If instead of the polygon ID you want to store in the mask the
thematic class to which the polygon belongs, you just add attr=iShp-
>IDLffShape::GetAttributes(i) in the beginning of the loop, replace in
the last line of the loop with mask[featis]= attr.(j) (where j is the
position in the dbf table of the field that contains the polygon
class).

Cheers

Guillermo

Guillermo writes: 

> ; populate the mask (assumes is in the same projection
> ;and covers the same extent as the geotiff)
> myshape->IDLffShape::GetProperty, N_ENTITIES=n
> FOR i=0L, n-1 DO BEGIN
>   feati=  myshape->IDLffShape::GetEntity(i)
>   featix= Round((Reform((*feati.vertices)[0,*])-x0)/psz)
>   featiy= Round((y0-Reform((*feati.vertices)[1,*]))/psz)
>   featis= POLYFILLV(featix, featiy, ns, nl)
>   IF featis[0] NE -1 THEN mask[featis]= feati.ishape +1
> ENDFOR

This is nearly identical to the solution I sent Paul
earlier today. The problem with IDLanROI is that it
can't be used in the native projected meter coordinates
that the image and shape files are in. The projected
meter coordinates have to be "converted" to "pixel"
coordinates by subtracting the offset and dividing
by the image range, before they can be loaded in the
object.

I didn't round my values, and I don't think you need
to do so here, even with PolyFillV. In fact, I think
you might get slightly more accurate values by not 
rounding, although this is a quibble with Paul's 
image.

I've made myself a note to write an article when
I get some time. :-)

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.")