Plot only on continent in matplotlib
I am drawing a map using basemap from matplotlib. The data are spreaded all over the world, but I just want to retain all the data on the continent and drop those on the ocean. Is there a way that I can filter the data, or is there a way to draw the ocean again to cover the data?
There's method in matplotlib.basemap: is_land(xpt, ypt)
It returns True
if the given x,y point (in projection coordinates) is over land, False
otherwise. The definition of land is based upon the GSHHS coastline polygons associated with the class instance. Points over lakes inside land regions are not counted as land points.
For more information, see here.
is_land()
will loop all the polygons to check whether it's land or not. For large data size, it's very slow. You can use points_inside_poly()
from matplotlib to check an array of points quickly. Here is the code. It doesn't check lakepolygons
, if you want remove points in lakes, you can add your self.
It took 2.7 seconds to check 100000 points on my PC. If you want more speed, you can convert the polygons into a bitmap, but it's a little difficult to do this. Please tell me if the following code is not fast enought for your dataset.
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.nxutils as nx
def points_in_polys(points, polys):
result = []
for poly in polys:
mask = nx.points_inside_poly(points, poly)
result.extend(points[mask])
points = points[~mask]
return np.array(result)
points = np.random.randint(0, 90, size=(100000, 2))
m = Basemap(projection='moll',lon_0=0,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
x, y = m(points[:,0], points[:,1])
loc = np.c_[x, y]
polys = [p.boundary for p in m.landpolygons]
land_loc = points_in_polys(loc, polys)
m.plot(land_loc[:, 0], land_loc[:, 1],'ro')
plt.show()
The HYRY's answer won't work on new versions of matplotlib (nxutils is deprecated). I've made a new version that works:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.path import Path
import numpy as np
map = Basemap(projection='cyl', resolution='c')
lons = [0., 0., 16., 76.]
lats = [0., 41., 19., 51.]
x, y = map(lons, lats)
locations = np.c_[x, y]
polygons = [Path(p.boundary) for p in map.landpolygons]
result = np.zeros(len(locations), dtype=bool)
for polygon in polygons:
result += np.array(polygon.contains_points(locations))
print result