Point in polygon test is one of the most useful functions when processing Geospatial data. Unfortunately, there is no official built-in function from Mathematica yet. Here is an example of showing usefulness of this function.
We have some data on Greenland icesheet thickness in format of {lon,lat, thickness}.
Our goal is to make a map to show the thickness of icesheet.
First try with ListDensityPlot:
You can see it is not working well, we need to limit the plot region inside the boundary of Greenland. Here is the place we can use the point-in-polygon test. In this example, I use inPolyQ2 from the answer by Simon Woods:
Still not right? What's wrong? The trick is to increase MaxPlotPoints to 100 at least:
ListDensityPlot[data, PlotRange -> All,
ColorFunction -> (ColorData["Rainbow"][1 - #] &),
RegionFunction -> (inPolyQ2[greenland_ploygon[[1, 1]], #1, #2] &),
MaxPlotPoints -> 150, MeshFunctions -> {#3 &}, Mesh -> 10]
It looks like a real map now.
2 comments:
This is very easy now with Mathematica V10:
In[1]:= poly = CountryData["Germany", "Polygon"];
In[2]:= coord = CityData["Berlin", "Coordinates"]
Out[2]= {52.52, 13.38}
In[3]:= RegionMember[poly /. GeoPosition -> Sequence, coord]
Out[3]= True
Nnevertheless, thanks for your sharing!
Kind regards, Armin
its very good post about Dentist .
SoftwareCorner
Post a Comment