define ROI from shapefile

[Brian McNoldy]

I am trying to calculate quantities averaged in US counties. I have a shapefile with the counties, but have not had much luck actually using that to do calculations (e.g., average some 2D data in each county in Colorado).  This seems like it should be fairly easy, but I'm coming up empty.  I assumed it would be related to a "region of interest" (ROI), but perhaps that was not correct.

Thanks for any tips!

Brian

[markus.sc...@gmail.com]

Did you try something like this? It should work, although you might have to add some keywords to computemask. Good Luck, Markus

roi =obj_new('IDLanROI',contour_line)

mask = roi.computemask() ne 0

average=mean(map_data[where(mask)])

[Brian McNoldy]

Thank you Markus,

I'm not quite sure what you mean; where does contour_line come from?  How would I pick and choose any or all counties in a state or even the whole country?  The region should be the whole county as an object, even if it has islands... would that still work? Sorry if I'm missing something obvious!

Brian

[markus.sc...@gmail.com]

In that case you have to isolate the border line for each island. contour_line has the format 2 x n. Then you can add or substract the masks to get enclaves and exclaves right. It's not obvious, but if you program and it step by step on real data it should be feasible.

[helder.m...@gmail.com]

Hi,

sorry if I'm a bit late...

I've played a bit around with shapefiles in the past. They are not my thing, but I had some fun with them.

Here you can see a demo I've made for european countries:

http://idl.marchetto.de/maps-and-european-countries/

Gettings statistics out of the boundary lines is pretty easy. Have a look at what David has done some time ago. He made a function called find_boundary():

http://www.idlcoyote.com/ip_tips/boundary.html

I hope it helps.

Cheers,

Helder

PS: I'm slowly abandoning IDL ... I had a good time with IDL. 20 years of fun, pain and successes... Somewhen I'll have to shutdown my not so much visited website ;-)

[Chang Liao]

Hi,

I wrote some code many moons ago to extract raster using shapefile, which might be relevant;

;;Read shapefile

oshp = OBJ_NEW('IDLffshape', sFilename_shapefile_in)

oshp -> GetProperty, n_entities = n_ent, Attribute_info = attr_info, $

n_attributes = n_attr, Entity_type = ent_type

roi_shp = LONARR(n_ent)

FOR ishp = 0, n_ent - 1 DO BEGIN

entitie = oshp -> GetEntity(ishp)

;;Check polygon

IF entitie.SHAPE_TYPE EQ 5 THEN BEGIN

record = *(entitie.VERTICES)

;;Convert coordinates

ENVI_CONVERT_FILE_COORDINATES, fid_in, xmap, ymap, record[0, *], record[1, *]

;;Create ROI

roi_shp[ishp] = ENVI_CREATE_ROI(ns = ns_in, nl = nl_in)

ENVI_DEFINE_ROI, roi_shp[ishp], /polygon, xpts = REFORM(xmap), ypts = REFORM(ymap)  

endfor

...

I wasn't using some latest APIs but they should be supported at the time of writing.

Good luck!

[markus.sc...@gmail.com]

Brian, I was not of  IDLffShape, but from reading on Helder's page, it seams to connect to what I wrote above about making a mask in the following way:

roi =obj_new('IDLanROI',pts[*,parts[k]:parts[k+1]-1])

mask = roi.computemask() ne 0

... combine masks to deal with islands ...

average=mean(map_data[where(mask)])

There might very well be a better way, but I'm not aware of it, as I'm only familiar with IDL and not with ENVI. 

Good Luck with the programming, Markus