nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

join tables based on location

ron_horne

Hi All

Can JMP join (update) tables based on location?

I have a data set of points with coordinates and a shapefile of polygons. I would like to update the dots data table with information from the shapefile based on location. So if a point is in a particular polygon it would update the table with the attributes of the polygon.

Thanks!

Ron

Accepted answers

ian_jmp

Coincidentally, I just found time to take a quick look. I make no claims for improved efficiency, but only offer it for comparison purposes. I did no checking or testing. For the subset you are using, it runs almost instantly. Note the OS/X paths.

NamesDefaultToHere(1);

ClearLog();

// Open dots file - this includes the observation we want to update

dots = Open( "$SAMPLE_DATA/Hurricanes.jmp" );

dots << delete rows( dots << get rows where( :Name( "Name and ID" ) != "NOTNAMED_675" ) );

// Open polygon source

polygons = Open( "/Library/Application Support/JMP/12/Maps/US-State-XY.jmp" );

// Open data source - this one includes the information to update

names = Open( "/Library/Application Support/JMP/12/Maps/US-State-Name.jmp" );

// Given a 'polygon matrix' pMat, and a shapeID and PartID therein, determines if the point (x, y) lies in the polygon defined by this shape and part

inShapePart =

Function({pMat, shapeID, partID, x, y}, {Default Local},

r = Loc(pMat[0,1] == shapeID & pMat[0,2] == partID);

InPolygon(x, y, pMat[r,3], pMat[r,4]);

);

// Given a 'polygon matrix' pMat, and a shapeID and PartID therein, determines if the point (x, y) lies in the rectangle inscribed by this shape and part

inShapePartRect =

Function({pMat, shapeID, partID, x, y}, {Default Local},

r = Loc(pMat[0,1] == shapeID & pMat[0,2] == partID);

maxX = Maximum(pMat[r, 3]);

minX = Minimum(pMat[r, 3]);

maxY = Maximum(pMat[r, 4]);

minY = Minimum(pMat[r, 4]);

if( (x <= maxX) & (x >= minX) & (y <= maxY) & (y >= minY), 1, 0);

);

pMat = polygons << getAsMatrix; // 'Polygon matrix'

x = Column(dots, "Longitude") << getValues; // x values

y = Column(dots, "Latitude") << getValues; // y values

loopLimits = polygons << Summary(Group( :Shape ID, :Part ID ), Invisible);

ns = Maximum(Column(loopLimits, "Shape ID") << getValues); // Number of shapes

np = Column(loopLimits, "Part ID") << getValues; // Number of parts in each shape

Close(loopLimits, NoSave);

For(s=1, s <= ns, s++, // Loop over shapes

For(p=1, p<=np[s], p++, // Loop over parts

For(d=NRow(x), d>=1, d--, // Loop over dots

// Is the dot inside the square that inscribes the polygon? (Relatively cheap!)

if (inShapePartRect(pMat, s, p, x[d], y[d]),

// If so, see if it's actually in the polygon (Relatively expensive!)

if (inShapePart(pMat, s, p, x[d], y[d]),

shapeName = Column(names, "Name")[s];

Print("Point ("||Char(x[d])||", "||Char(y[d])||") found in Shape "||Char(s)||", Part "||Char(p)||", with Name "||shapeName||".");

// NB: As each dot is 'found', delete the corresponding entries in 'x' and'y' to try to speed things up

x[d] = []; y[d] = []

);

Close(dots, NoSave);

Close(polygons, NoSave);

Close(names, NoSave);

All comments

ian_jmp

I think you will need scripting. The 'InPolygon()' JSL function should be useful.

ron_horne

Thanks Ian@JMP

This looks promising. will go through what i need it to do and see if it works.

ron

msharp

I would love to see a sample code if you get it to work.

ron_horne

Hi All,

I have manage to get the script working. It takes each point and matches it with the polygon it is surrounded by to join the tables. In this case it adds the name of the state to each dot.

The script does its job but perhaps not in the most efficient way. Anny ideas for improved efficiency are welcome.

If you get to use it please share your experience.

Clear Globals();

Names Default To Here( 1 );

// open dots file - this includes the observation we want to update

dots = Open( "C:\Program Files (x86)\SAS\JMPSW\11\Samples\Data\Hurricanes.jmp" );

// reduce the size of the data for the illustration

dots << delete rows( dots << get rows where( :Name( "Name and ID" ) != "NOTNAMED_675" ) );

// open polygon source

polygons = Open( "C:\Program Files (x86)\SAS\JMPSW\11\Maps\US-State-XY.jmp", Private );

// open data source - this one includes the information to update

data = Open( "C:\Program Files (x86)\SAS\JMPSW\11\Maps\US-State-Name.jmp", Private );

// introduce the column with the attributes to be updated in the dots table

dots << New Column( "State", Character, Nominal );

dots << Set Label Columns( :State );

Wait( 1 );

// extract the polygons into matrices to improve efficiency

XMATRIX = polygons:X << get as matrix; // x values of polygons

YMATRIX = polygons:Y << get as matrix; // y value of polygons

Start = Tick Seconds();

// run through the dots across the polygons and check for matches.

For( i = 1, i <= N Rows( dots ), i++, // for each dot

For( ii = 1, ii <= Col Max( polygons:Shape ID ), ii++, // for each shape ID

ROWLISTII = polygons << get rows where( polygons:Shape ID == ii ); // get the rows for each Shape ID

For( iii = 1, iii <= Max( Polygons:Part ID[ROWLISTII] ), iii++, // for each Part ID according to the number of Part IDs for each Shape ID

// get the rows for each polygon

rowsii = polygons << get rows where( And( polygons:Shape ID == ii, polygons:Part ID == iii ) );

// get the X values for each polygon

XMATRIXROWSII = XMATRIX[rowsii];

// get the Y values for each polygon

YMATRIXROWSII = YMATRIX[rowsii];

// compare the location of the dot with the polygon

result = In Polygon( dots:longitude[i], dots:Latitude[i], XMATRIXROWSII, YMATRIXROWSII );

If( result == 1, Show( i, ii, data:Name[ii] ) ); // this will send the list of matching pairs to the log

If( result == 1,

dots:State[i] = data:Name[ii];// this will make the update of the attribute from the polygon information (data) to the dot (dots)

Break(); // this will make the third loop (iii) stop checking more polygons of the same Shape ID.

);

If( result == 1, Break() ); // this will make the second loop (ii) stop checking for more Shape ID's for the same dot and move to the next dot in the list.

)

);

End = Tick Seconds();

Runtime = End - Start;

Show( runtime );

// Wait( 1 );

dots << Bivariate(

Y( :Latitude ),

X( :Longitude ),

SendToReport(

Dispatch( {}, "1", ScaleBox, {Scale( "Geodesic" ), Min( -94.0515072337543 ), Max( -60.9484927662457 ), Inc( 5 ), Minor Ticks( 1 )} ),

Dispatch( {}, "2", ScaleBox, {Scale( "Geodesic" ), Min( 7.96034811422145 ), Max( 49.6893151267156 ), Inc( 5 ), Minor Ticks( 1 )} ),

Dispatch(

{},

"Bivar Plot",

FrameBox,

{Frame Size( 449, 566 ), Background Map( Boundaries( "US States" ) ), Grid Line Order( 3 ), Reference Line Order( 4 ),

Row Legend(

State,

Color( 1 ),

Color Theme( "JMP Default" ),

Marker( 0 ),

Marker Theme( "" ),

Continuous Scale( 0 ),

Reverse Scale( 0 ),

Excluded Rows( 0 )

)}

)

);

Map1.png

ian_jmp

NamesDefaultToHere(1);

ClearLog();

// Open dots file - this includes the observation we want to update

dots = Open( "$SAMPLE_DATA/Hurricanes.jmp" );

dots << delete rows( dots << get rows where( :Name( "Name and ID" ) != "NOTNAMED_675" ) );

// Open polygon source

polygons = Open( "/Library/Application Support/JMP/12/Maps/US-State-XY.jmp" );

// Open data source - this one includes the information to update

names = Open( "/Library/Application Support/JMP/12/Maps/US-State-Name.jmp" );

// Given a 'polygon matrix' pMat, and a shapeID and PartID therein, determines if the point (x, y) lies in the polygon defined by this shape and part

inShapePart =

Function({pMat, shapeID, partID, x, y}, {Default Local},

r = Loc(pMat[0,1] == shapeID & pMat[0,2] == partID);

InPolygon(x, y, pMat[r,3], pMat[r,4]);

);

// Given a 'polygon matrix' pMat, and a shapeID and PartID therein, determines if the point (x, y) lies in the rectangle inscribed by this shape and part

inShapePartRect =

Function({pMat, shapeID, partID, x, y}, {Default Local},

r = Loc(pMat[0,1] == shapeID & pMat[0,2] == partID);

maxX = Maximum(pMat[r, 3]);

minX = Minimum(pMat[r, 3]);

maxY = Maximum(pMat[r, 4]);

minY = Minimum(pMat[r, 4]);

if( (x <= maxX) & (x >= minX) & (y <= maxY) & (y >= minY), 1, 0);

);

pMat = polygons << getAsMatrix; // 'Polygon matrix'

x = Column(dots, "Longitude") << getValues; // x values

y = Column(dots, "Latitude") << getValues; // y values

loopLimits = polygons << Summary(Group( :Shape ID, :Part ID ), Invisible);

ns = Maximum(Column(loopLimits, "Shape ID") << getValues); // Number of shapes

np = Column(loopLimits, "Part ID") << getValues; // Number of parts in each shape

Close(loopLimits, NoSave);

For(s=1, s <= ns, s++, // Loop over shapes

For(p=1, p<=np[s], p++, // Loop over parts

For(d=NRow(x), d>=1, d--, // Loop over dots

// Is the dot inside the square that inscribes the polygon? (Relatively cheap!)

if (inShapePartRect(pMat, s, p, x[d], y[d]),

// If so, see if it's actually in the polygon (Relatively expensive!)

if (inShapePart(pMat, s, p, x[d], y[d]),

shapeName = Column(names, "Name")[s];

Print("Point ("||Char(x[d])||", "||Char(y[d])||") found in Shape "||Char(s)||", Part "||Char(p)||", with Name "||shapeName||".");

// NB: As each dot is 'found', delete the corresponding entries in 'x' and'y' to try to speed things up

x[d] = []; y[d] = []

);

Close(dots, NoSave);

Close(polygons, NoSave);

Close(names, NoSave);

ron_horne

8 and a half hours later the whole file of hurricanes was joined. since most of the dots are in the sea the script took some time to verify that each one of them was not in any of the 52 states (84 polygons).

ron_horne

thank you Ian@JMP,

this is way better than my approach and much faster.

while producing the loops in my script i noticed that the XY table is not perfectly clean and your approach is also robust against that.(shape ID 10 has part ID 1,3,4,5 - skipping 2 ).

ron

msharp

To further speed the script up, you can check if the point is in "America" square first. Appending onto Ian's script

inAmericaRect =

Function({pMat, x, y}, {Default Local},

maxX = Maximum(pMat[0, 3]);

minX = Minimum(pMat[0, 3]);

maxY = Maximum(pMat[0, 4]);

minY = Minimum(pMat[0, 4]);

if( (x <= maxX) & (x >= minX) & (y <= maxY) & (y >= minY), 1, 0);

);

...

// Is the dot inside the Square of America?

if (inAmericaRect(pMat, x, y),

Of course, given Puerto Rico, this could potentially slow down the script by adding an extra if statement. I still feel there is enough data outside to warrant it. But to fix that, you could simply split out Puerto Rico or other islands and add an addition inPuertoRicoRect function. You could do additional filtering for other countries using world-name.jmp and world-xy.jmp. Or you could use the coastline-Name.jmp and coastline-xy.jmp to first determine if it's on land. But, I'm sure at that point your script will be going outside the scope of the project.