is_land slow for multiple coordinates

[jakobloeg]

Hello, I have an application where I would like to call is_land multiple times and for many coordinates, in the order of millions. Is there a way to speed up this process, instead of having to call is_land with one coordinate at the time. This is simply too slow for my application

[WeatherGod]

I doubt it. Not without taking advantage of some spatial indexing code such that geopandas or (I think) shapely can take advantage of if explicitly set. Btw, in case you aren't aware, basemap is considered deprecated, and you should move off of it, moving towards packages like cartopy.

…

On Wed, May 29, 2019 at 2:33 PM jakobloeg ***@***.***> wrote: Hello, I have an application where I would like to call is_land multiple times and for many coordinates, in the order of millions. Is there a way to speed up this process, instead of having to call is_land with one coordinate at the time. This is simply too slow for my application — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#462?email_source=notifications&email_token=AACHF6FNHLKA47Y7PQD6I5DPX3D7JA5CNFSM4HQQKKNKYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4GWR2QWA>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AACHF6GYQJH6UWHV2UKZT6TPX3D7JANCNFSM4HQQKKNA> .

[molinav]

You can create a vectorized version of is_land with the help of maskoceans:

def is_land(lon, lat, resolution="l", grid=5):

    from mpl_toolkits.basemap import maskoceans

    lon, lat = np.broadcast_arrays(lon, lat)
    dummy = np.zeros(lon.shape, dtype=bool)
    mask = ~maskoceans(lon, lat, dummy, resolution=resolution, grid=grid).mask

    return mask

and use it e.g. as follows:

step = 0.25
lons = +0.5 * step + np.arange(-180, +180, +step)[None, :]
lats = -0.5 * step + np.arange(+90, -90, -step)[:, None]

land = is_land(lons, lats)

I use a similar function for my own purposes, and it is way faster than the normal is_land. However, the maskoceans function from my current version of basemap (1.0.8) cannot return correct values in the neighbourhood of latitudes equal to ±90° or longitudes equal to ±180°. This seems to happen when you request a coordinate that is outside the limits from the underlying GSHHG maps. However, Basemap.is_land is capable of handling such cases. I should report it at some point.