X
PrevPrev Go to previous topic
NextNext Go to next topic
Last Post 14 Mar 2022 11:00 AM by  Ben Castellani
How to determine which US state a point is located in?
 0 Replies
Sort:
You are not authorized to post a reply.
Author Messages

Ben Castellani



Basic Member


Posts:130
Basic Member


--
14 Mar 2022 11:00 AM
    There is no routine built into IDL to do this by default. However, it is very simple using the stock states shapefile and the IDLanROI object class: https://www.l3harrisgeospatial.com/docs/idlanroi.html

    See the provided example function below which can be used to identify the state containing the specified point:

    function stateFinder, longitude, latitude

    ; NAME:
    ; stateFinder
    ;
    ; PURPOSE:
    ; This function returns the name of the US state containing the specified point.
    ;
    ; CALLING SEQUENCE:
    ; Result = stateFinder(longitude,latitude)
    ;
    ; INPUTS:
    ; LONGITUDE: A scalar containing the longitude of the point.
    ; LATITUDE: A scalar containing the latitude of the point.
    ;
    ; RETURN VALUE:
    ; A string containing the name of the US state containing the point. If the point lies outside all states, 'None' is returned.

    ;Load the states shapefile from the IDL installation
    state_shapefile= FILEPATH('states_high.shp', SUBDIRECTORY=['resource','maps','shape'])
    myshape1 = OBJ_NEW('IDLffShape',state_shapefile)

    ;Get the number of polygons in file
    myshape1->GetProperty, N_ENTITIES=num_ent

    FOR x=0L, (num_ent-1) DO BEGIN

    ; Get the Attributes for the active state
    attr = myshape1->GetAttributes(x)

    ;Store state name
    state_name=attr.ATTRIBUTE_0

    ; Read the vertices of the state
    ent = myshape1->GetEntity(x)
    xpoints = (*ent.VERTICES)[0, *]
    ypoints = (*ent.VERTICES)[1, *]

    ;Determine whether the specified point is within this state
    oROI = OBJ_NEW( 'IDLanROI',xpoints,ypoints)
    testInside = oROI -> ContainsPoints(longitude,latitude)

    ;If state is found, stop searching states
    if testInside gt 0 then break
    ENDFOR

    ;If no states contain the point, return 'None'
    if testInside eq 0 then state_name='None'
    RETURN,state_name

    END


    Here are some example outputs:

    IDL> stateFinder(-105,40)
    Colorado
    IDL> stateFinder(-75,40)
    New Jersey
    IDL> stateFinder(-90,25)
    None
    You are not authorized to post a reply.