[GIS] Masking land with matplotlib basemap in Python

geolocationmatplotlib-basemappython

I've been trying to mask land with .is_land from mpl_toolkits.basemap.

When running following code:

map = Basemap(llcrnrlon = 7.298914095230308, llcrnrlat = 58.98235690451632, 
           urcrnrlon = 12.27072348324015, urcrnrlat = 
           57.92306182768044,projection='cyl', resolution='f')
value = map.is_land(11.61168822665539, 57.86868795294363)

print(value) returns False

but when not specifying region:

map = Basemap(projection='cyl', resolution='f')
value = map.is_land(11.61168822665539, 57.86868795294363)

print(value) returns True

I can't figure out why.

I really do need to specify region or else the code will be 5-6 times slower to run.

Best Answer

Basemap.is_land() uses point-in-polygon algorithm to find a solution, it requires one of the GSHHS coastline polygons in the data files supplied with Basemap package. With resolution='f', you ask for the biggest file (80 MB) with most detailed coastlines and lakes. When you don't specify lat-long extents, the whole world coastline data is loaded and run. But, when lat-long limits are specified, the world coastline polygons are sliced/cut for the required parts. This simplifies the coastline data before use, some vertex points might be deleted/added, causing deviation on parts of the coastlines/lakes. It also filters out polygons outside specified extents. This process could be the reasons to your findings. Also note that points over lakes inside land regions are not counted as land points.

Related Solutions

[GIS] Shapefile projection for Matplotlib Basemap

Use directly the GDAL module rather than using it through (Geo)Django:

from osgeo import osr
ori =  osr.SpatialReference()
desti =  osr.SpatialReference()
ori.ImportFromProj4("+init=EPSG:2227")
desti.ImportFromProj4("+init=EPSG:4326")
transformation = osr.CoordinateTransformation(ori,desti)
transformation1.TransformPoint(5979385.3163340045, 2085856.3243659574)
(-122.5129551683127, 37.706167841402952, 0.0)

or the pyproj module (with preserve_units=True because pyproj assumes that your coordinates are in meters and this is not the case for EPSG:2227):

from pyproj import Proj, transform
ori = Proj(init='epsg:2227',preserve_units=True) 
dest= Proj(init='EPSG:4326',preserve_units=True)
transform(ori, dest,x,y)
(-122.51295516831273, 37.706167841402959)

Combine with Fiona:

from pyproj import Proj, transform
import fiona
from fiona.crs import from_epsg
with fiona.open("/data/sfpd_districts.shp") as shp
    ori = Proj(shp.crs,preserve_units=True ),
    dest= Proj(init='EPSG:4326',preserve_units=True)
    with fiona.open('/data/sfpd_districtsWGS84.shp', 'w', 'ESRI Shapefile', shp.schema.copy(), crs=from_epsg(4326) as output:
        for point in shp:
            x,y =  point['geometry']['coordinates']
            point['geometry']['coordinates'] = transform(ori, dest,x,y)
            output.write(point)

[GIS] How accurate is Pythons is_land

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

Best Answer

Related Solutions

[GIS] Shapefile projection for Matplotlib Basemap

[GIS] How accurate is Pythons is_land

Related Question