Need to detect map's land/ocean points with higher accuracy than maskoceans(grid=1.25,resolution='f')

Question

I have a grid of points on the map and I want to detect which point is in land or ocean. I have tried to do this in 3 ways but non are satisfactory.

1. Using maskoceans, a simple example is:

   from mpl_toolkits.basemap import maskoceans
   import numpy
   # my grid is given by lons and lats 1D arrays.
   Array = numpy.ma.zeros((lats.size,lons,size))
   Lons,Lats = numpy.meshgrid(lons,lats)
   MaskedArray = maskoceans(Lons,Lats,Array,resolution='f',grid=1.25)
   

Now I have MaskedArray where locations on ocean are masked. This method is really fast compared to others, but inaccurate near the coast for small scales. In the following plot, red dots are identified as land, black points are identified as ocean. As you see even if I have used the resolution='f' (the full reslution) and smallest masked grid (grid=1.25) there are erounous points in ocean that are marked as land in red. It seems that the f resolution of maskoceans is not in the same order of the f resolution of Basemap.

Image of points detected with maskoceans

2. Use basemap.is_land(). Example is:

   # Create a basemap object, this takes a while!
   MinLat = numpy.min(lats) - 0.01
   MinLon = numpy.min(lons) - 0.01
   MaxLat = numpy.max(lats) + 0.01
   MaxLon = numpy.max(lons) + 0.01
   MidLat = numpy.mean(lats)
   MidLon = numpy.mean(lons)
   from mpl_tookits.basemap import Basemap
   map = Basemap(projection='aeqd',llcrnrlat=MinLat,llcrnrlon=MinLon,urcrnrlat=MaxLat,urcrnrlon=MaxLon,area_thresh=0.01,lon_0=MidLon,lat_0=MidLat,resolution='f')
   Array = numpy.zeros((lats.size,lons.size),dtype=bool)
   for j in range(lats.size):
       for i in range(lons.size):
           x,y = map(lon[i],lat[j])
           if map.is_land(x,y):
               Array[j,i] = True
   

This is the most accurate method. As you see in the figure below, all points are detected correctly. The big issue with is_land() is that it is very slow. For large grids this method takes minutes (in order of an hour for very large grids), which maskoceans take less than a second.

Image of points detected with Basemap.is_land

3. Using Polygon.contains_point. See an example here. This does not work for me, since the polygons of the coastlines should be closed in order to use conians_point. Hence the map object should be defined on the the full globe with

    map = Basemap(projection='aeqd',lon_0=0,resolution='f')
   

If the map is defined on a window of lon and lat, the coastline polygon is not closed. But the loading the full map with the full resolution takes a very long time and memory, and detecting even one point takes forever.

Is there an alternative method to identify land/ocean points on a given grid? It would prefer to use maskoceans if the resolution issue is fixed. It is possible to increase the resolution of maskoceans() up to the resolution of Basemap?