Determine if map coordinates are within area

[EvonS]

I have this information below that define areas on a map. Is there a way to determine if a longitude and latitude are within one of these areas?

OBJECTID	AGENCY	DIV	UNIT_NAME	ADDRESS	CITY	AREA_SQKM	Shape__Area	Shape__Length
1​	Toronto Police Service	D11	11 Division	2054 Davenport Rd	Toronto	18.67845803​	35692935.07​	28190.11106​
2​	Toronto Police Service	D12	12 Division	200 Trethewey Dr	Toronto	24.78034332​	47417573.68​	33029.13457​
3​	Toronto Police Service	D13	13 Division	1435 Eglinton Av W	Toronto	18.60191063​	35591722.93​	24616.57917​
4​	Toronto Police Service	D14	14 Division	350 Dovercourt Rd	Toronto	14.17525415​	27083979.14​	30017.70564​
5​	Toronto Police Service	D22	22 Division	3699 Bloor St W	Toronto	68.12188386​	130097609.1​	72765.03189​
6​	Toronto Police Service	D23	23 Division	5230 Finch Av W	Toronto	57.64228008​	110369052.3​	54270.05361​
7​	Toronto Police Service	D31	31 Division	40 Norfinch Dr	Toronto	42.66749447​	81789149.85​	44477.14656​
8​	Toronto Police Service	D32	32 Division	30 Ellerslie Av	Toronto	61.06048027​	117079222.9​	43658.49725​
9​	Toronto Police Service	D33	33 Division	50 Upjohn Rd	Toronto	52.7888815​	101244281.9​	45148.84874​
10​	Toronto Police Service	D41	41 Division	2222 Eglinton Av E	Toronto	44.57511239​	85395988.8​	58875.05415​
11​	Toronto Police Service	D42	42 Division	242 Milner Av E	Toronto	84.61174885​	162528804.9​	59824.09155​
12​	Toronto Police Service	D43	43 Division	4331 Lawrence Av E	Toronto	59.37000474​	113873893.7​	54424.50711​
13​	Toronto Police Service	D51	51 Division	51 Parliament St	Toronto	8.601682841​	16437802.1​	25324.5062​
14​	Toronto Police Service	D52	52 Division	255 Dundas St W	Toronto	9.228055772​	17625448.21​	35966.54447​
15​	Toronto Police Service	D53	53 Division	75 Eglinton Av W	Toronto	31.29485302​	59895251.81​	37395.46015​
16​	Toronto Police Service	D55	55 Division	101 Coxwell Av	Toronto	41.03398633​	78486159.8​	82307.80296​

 

[EvonS]

@KRice and @lrobbo314 , thank you so much. I've incorporated your work into my submission form and it works perfectly. And thank you to everyone else who took the time to reply. You all are awesome.

 

[KRice]

That's good news. We're happy to help.

 

[lrobbo314]

@KRice,

Nice workbook you uploaded. I found the same geojson files and played around with the PowerQuery portion of it.

I managed to get it down to just the 2 queries in PQ.

Here is what I came up with.

It appends a letter to Divisions that had multiple polygons.

This was a really fun project and learned some cool new things along the way.

tbl_Main

let
    Source = Json.Document(File.Contents("C:\Users\UserName\Downloads\TPS_Police_Divisions.geojson")),
    Custom1 = Source[features],
    ToTable = Table.FromList(Custom1, Splitter.SplitByNothing(), null, null, ExtraValues.Error),
    Main = Table.ExpandRecordColumn(ToTable, "Column1", {"type", "properties", "geometry"}, {"type", "properties", "geometry"}),
    Properties = Table.ExpandRecordColumn(Main, "properties", {"OBJECTID", "AGENCY", "DIV", "UNIT_NAME", "ADDRESS", "CITY", "AREA_SQKM", "Shape__Area", "Shape__Length"}, {"OBJECTID", "AGENCY", "DIV", "UNIT_NAME", "ADDRESS", "CITY", "AREA_SQKM", "Shape__Area", "Shape__Length"}),
    Geometry = Table.ExpandRecordColumn(Properties, "geometry", {"type", "coordinates"}, {"type.1", "coordinates"}),
    Coordinates = Table.ExpandListColumn(Geometry, "coordinates"),
    Index = Table.AddIndexColumn(Coordinates, "Index", 1, 1, Int64.Type),
    NewDiv = Table.AddColumn(Index, "DIVX", (ro)=> if List.Count(List.Select(Index[DIV], each _ =  ro[DIV])) = 1 then ro[DIV] else ro[DIV] & "-" & Character.FromNumber(List.Count(List.Select(List.FirstN(Index[DIV],ro[Index]), each _ = ro[DIV]))+64)),
    RC = Table.RemoveColumns(NewDiv,{"type", "OBJECTID", "DIV", "coordinates", "Index", "type.1"}),
    Reorder = Table.ReorderColumns(RC,{"AGENCY", "DIVX", "UNIT_NAME", "ADDRESS", "CITY", "AREA_SQKM", "Shape__Area", "Shape__Length"}),
    Types = Table.TransformColumnTypes(Reorder,{{"Shape__Length", type number}, {"Shape__Area", type number}, {"AREA_SQKM", type number}, {"CITY", type text}, {"ADDRESS", type text}, {"UNIT_NAME", type text}, {"DIVX", type text}, {"AGENCY", type text}})
in
    Types

tbl_Coordinates

let
    Source = Json.Document(File.Contents("C:\Users\UserName\Downloads\TPS_Police_Divisions.geojson")),
    Custom1 = Source[features],
    ToTable = Table.FromList(Custom1, Splitter.SplitByNothing(), null, null, ExtraValues.Error),
    Main = Table.ExpandRecordColumn(ToTable, "Column1", {"type", "properties", "geometry"}, {"type", "properties", "geometry"}),
    Properties = Table.ExpandRecordColumn(Main, "properties", {"OBJECTID", "AGENCY", "DIV", "UNIT_NAME", "ADDRESS", "CITY", "AREA_SQKM", "Shape__Area", "Shape__Length"}, {"OBJECTID", "AGENCY", "DIV", "UNIT_NAME", "ADDRESS", "CITY", "AREA_SQKM", "Shape__Area", "Shape__Length"}),
    Geometry = Table.ExpandRecordColumn(Properties, "geometry", {"type", "coordinates"}, {"type.1", "coordinates"}),
    Coordinates = Table.ExpandListColumn(Geometry, "coordinates"),
    Index = Table.AddIndexColumn(Coordinates, "Index", 1, 1, Int64.Type),
    NewDiv = Table.AddColumn(Index, "DIVX", (ro)=> if List.Count(List.Select(Index[DIV], each _ =  ro[DIV])) = 1 then ro[DIV] else ro[DIV] & "-" & Character.FromNumber(List.Count(List.Select(List.FirstN(Index[DIV],ro[Index]), each _ = ro[DIV]))+64)),
    ExpandMult = Table.TransformColumns(NewDiv,{{"coordinates", each if List.Count(_) = 1 then _{0} else _}}),
    RC = Table.RemoveColumns(ExpandMult,{"type", "AGENCY", "ADDRESS", "CITY", "AREA_SQKM", "Shape__Area", "Shape__Length", "type.1", "OBJECTID", "Index", "DIV"}),
    Reorder = Table.ReorderColumns(RC,{"DIVX", "UNIT_NAME", "coordinates"}),
    ExpandedCoords = Table.ExpandListColumn(Reorder, "coordinates"),
    ExtractValues = Table.TransformColumns(ExpandedCoords, {"coordinates", each Text.Combine(List.Transform(_, Text.From), ";"), type text}),
    Split = Table.SplitColumn(ExtractValues, "coordinates", Splitter.SplitTextByDelimiter(";", QuoteStyle.Csv), {"Longitude", "Latitude"}),
    Types = Table.TransformColumnTypes(Split,{{"Latitude", type number}, {"Longitude", type number}, {"UNIT_NAME", type text}, {"DIVX", type text}})
in
    Types

CAN InPoly II.xlsm
	A	B	C	D	E	F	G	H	I	J
1					Bounding Boxes
2		Street Address	City		Division	LongMin	LongMax	LatMin	LatMax	In Polygon
3	User Inputs	45 Esandar Dr	Toronto		D11	-79.5122	-79.4286	43.63171	43.67511	Not In Polygon
4	Lat-Long of address	43.706916	-79.356874		D12	-79.5396	-79.456	43.66264	43.72469	Not In Polygon
5	Location in:	D53			D13	-79.4698	-79.4073	43.66727	43.71927	Not In Polygon
6	Send to:	75 Eglinton Av W	Toronto		D14	-79.4462	-79.3916	43.62627	43.67544	Not In Polygon
7					D22	-79.6087	-79.4691	43.581	43.68413	Not In Polygon
8					D23-A	-79.6393	-79.5102	43.66463	43.76303	Not In Polygon
9					D23-B	-79.5654	-79.5654	43.73736	43.73743	Not In Polygon
10					D31	-79.5852	-79.4725	43.71282	43.7839	Not In Polygon
11					D32	-79.4845	-79.3811	43.71081	43.80327	Not In Polygon
12					D33	-79.3963	-79.3026	43.71624	43.81566	Not In Polygon
13					D41	-79.3198	-79.2191	43.67107	43.77934	Not In Polygon
14					D42	-79.3413	-79.1339	43.76837	43.85546	Not In Polygon
15					D43	-79.2648	-79.1155	43.70829	43.80463	Not In Polygon
16					D51	-79.3868	-79.3377	43.63321	43.67533	Not In Polygon
17					D52-A	-79.4049	-79.3494	43.61196	43.63535	Not In Polygon
18					D52-B	-79.4039	-79.3742	43.63702	43.67023	Not In Polygon
19					D53	-79.4298	-79.3354	43.66667	43.73037	In Polygon
20					D55	-79.3664	-79.28	43.61274	43.7263	Not In Polygon
User_Input

Cell Formulas
Range	Formula
E3:E20	E3	=UNIQUE(Table_tbl_Main[DIVX])
F3:F20	F3	=MINIFS(Table_tbl_Coordinates[Longitude],Table_tbl_Coordinates[DIVX],User_Input!E3#)
G3:G20	G3	=MAXIFS(Table_tbl_Coordinates[Longitude],Table_tbl_Coordinates[DIVX],User_Input!E3#)
H3:H20	H3	=MINIFS(Table_tbl_Coordinates[Latitude],Table_tbl_Coordinates[DIVX],User_Input!E3#)
I3:I20	I3	=MAXIFS(Table_tbl_Coordinates[Latitude],Table_tbl_Coordinates[DIVX],User_Input!E3#)
B4:C4	B4	=TEXTSPLIT(TEXTBEFORE(TEXTAFTER(WEBSERVICE("http://geocoder.ca/?locate=" & SUBSTITUTE(B3," ","+")&"%2c+"&C3&"%2c+Ontario&geoit=GeoCode+it!"),"""ICBM"" content="""),""""),", ")+0
B5	B5	=INDEX(E3#,MATCH("In Polygon",$J$3:$J$20,0))
B6:C6	B6	=FILTER(Table_tbl_Main[[ADDRESS]:[CITY]],Table_tbl_Main[DIVX]=User_Input!B5)
J3:J20	J3	=IF(PtInPoly($C$4,$B$4,FILTER(Table_tbl_Coordinates[[Longitude]:[Latitude]],Table_tbl_Coordinates[DIVX]=User_Input!E3)),"","Not ")&"In Polygon"
Dynamic array formulas.

And yes, Rick's PtInPoly algorithm works perfectly.

 

[lrobbo314]

Also,

Here is a conversion of Rick's PtInPoly algorithm to a LAMBDA function.

=LAMBDA(XC,YC,P,
    IF(MOD(SUM(
        SCAN(0,SEQUENCE(ROWS(P)-1),
            LAMBDA(s,c,
                IF(XOR(INDEX(P,c,1)>XC,INDEX(P,c+1,1)>XC),
                    LET(
                        m,(INDEX(P,c+1,2)-INDEX(P,c,2))/(INDEX(P,c+1,1)-INDEX(P,c,1)),
                        b,(INDEX(P,c,2)*INDEX(P,c+1,1)-INDEX(P,c,1)*INDEX(P,c+1,2))/(INDEX(P,c+1,1)-INDEX(P,c,1)),
                        IF(m*XC+b>YC,s+1,s+0)
                    ),0
                )
            )
        )
    ),2),"In Polygon","Not In Polygon"
    )
)

 

[lrobbo314]

Just curious @Rick Rothstein ... Did you come up with that and know how to do the math off the top of your head? Or did you look up pseudo-code and translate it to VBA?

 

[Rick Rothstein]

Just curious @Rick Rothstein ... Did you come up with that and know how to do the math off the top of your head? Or did you look up pseudo-code and translate it to VBA?

I developed the original code for that function many years ago (maybe 15 to 20 years ago now) when I was volunteering in the Visual Basic (the original compiled version, not the one Microsoft calls Visual Basic now) newsgroups (the predecessor to forum for you young'un out there). I developed the math myself... I was a math major in college and a licensed Professional Civil Engineer (designing roads and such) for the bulk of my working career, so developing the function was not that difficult for me. If you look through the entire thread at the link I provided, you will find some modification on the code in response to requests that were made to me.

 

[lrobbo314]

Wow. That is very impressive. I tried to figure out what was going on there, and it's a bit beyond me. But I'm inspired to go dig deeper into it.

Thanks for the update.

 

[KRice]

@lrobbo314, I came across this reference when I was investigating various "point in polygon" strategies...before Rick pointed toward a very nice solution that could be incorporated into the OP's form.

Geometry Algorithms Home

At one point, that page offered the entire book, but the book proper appears to have been removed and is now available through other outlets (I have a pdf version of it). The book contains a couple of strategies for point-in-polygon and some good references, as well as strategies for solving many other types of geometry problems. While the book contains a lot of pseudo-code, the C++ code for the algorithms is still available for free download through the website. One of the strategies describes counting edge crossings.
@Rick Rothstein, you mentioned 15-20 years ago. After exploring the link you provided, I found a download for a file you created in 2003 that included the point-in-polygon UDF. It's remarkable that a generalized algorithm developed long ago has worked so well for solving a variety of different problems.

 

[EvonS]

I'm glad my post inspired so much great work. If anyone is curious, I've attached the file where I put your formulas to practical use. Let me know what you think. FOCUS Agency Submission Form 10.0.xlsm

 

[Rick Rothstein]

What was your search criteria for finding that 2003 article? The earliest I could find was one written 2013. What really amazed me across the years about the Point In Polygon function was the variety of areas people wrote that they found a use for it in... an eye doctor wrote that he used it for some kind of vision mapping need, someone used it to determine the spacing of streetlights within a section of their town, another used it for some kind of soil mapping application, someone even used it for determining the winner of a hot air balloon contest, and the list goes on.