How accurate is Pythons is_land?

Question

I was planning to use Pythons is_land function, to differentiate between land and sea. However, after some inspection, I have major concerns in the accuracy. I can't find anything about this in the documentation, so I hoped anyone here could help me.

I initiate the basemap with high resolution, and a low area threshhold to help differentiate:

map = Basemap(projection='merc', resolution = 'h', area_thresh=0.001)

When I try this function, it returns true even though I choose points clearly not on land:

>>> map.is_land(39.784004,3.345337)
True
>>> map.is_land(39.789808,3.3395)
True
>>> map.is_land(39.799041,3.332291)
True
>>> map.is_land(39.863371,3.346367)
True

The three last ones is out at sea.

Anyone here that has used this function, and maybe got some tips to how I can make it more accurate, or if there is some (preferably) Python substitution, or other open source programmable solution to check if a number of points is at land or not.

Thought maybe I needed to map the points first, no better:

>>> map = Basemap(projection='merc', resolution = 'h', area_thresh=0.001)
>>> lat,long = map(39.784004,3.345337)
>>> map.is_land(lat,long)
True
>>> lat,long = map(39.789808,3.3395)
>>> map.is_land(lat,long)
True
>>> lat,long = map(39.799041,3.332291)
>>> map.is_land(lat,long)
True
>>> lat,long = map(39.863371,3.346367)
>>> map.is_land(lat,long)
True

Answer 1

The reason why it doesn't work is because the arguments for is_land(xin, yin) are referencing the grid, and aren't latitudinal and longitudinal coordinates, they're simply the index values.

Since is_land(xin, yin) takes index values and not lat/long coordinates, you need to find out the index values on the map for the coordinates you're interested in.

map = basemap.Basemap(
  area_thresh=10,
  resolution="l",
  llcrnrlon=0.,
  llcrnrlat=-80.,
  urcrnrlon=360,
  urcrnrlat=80
)
lat, long = -37, 158
x,y = map(long, lat)
>>> map.is_land(x,y)
False

Those lat and long coordinates are off the coast of Australia, so this works. x and y don't need to be integers, they can be floats.

EDIT:

Notice that the arguments for map are map(long, lat).

Answer 2

The documentation says "Calling a Basemap class instance with the arguments lon, lat will convert lon/lat (in degrees) to x/y map projection coordinates (in meters)".

So you are passing in points around 39.789808E, 3.3395N - those are definitely on land (in Kenya). If you were hoping for points near the coast in Spain, you've just reversed the coordinates. Read the documentation again - "lon, lat".

You are using lat and long as variable names, but you are really getting X and Y (which would be closer to longitude and latitude, not latitude and longitude).

This is a very common issue in geospatial applications. It usually helps to remember its mostly Cartesian concepts, so X then Y.

Answer 3

I've written a small package with a function that is much faster than basemap.is_land() https://github.com/toddkarin/global-land-mask

It provides a global binary mask of land/ocean at 1 km resolution.