# How accurate is Pythons is_land?

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
``````

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)`.

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.

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.