[GIS] How accurate is Pythons is_land

basemappython

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

Best Answer

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

Related Question