Check if Lat/Lon point is over land or ocean

0
votes

I have an array of latitudes (lats, min=-88, max=88, shape=89) and longitudes (lons, min=0, max=358, shape=180), and a land mask (land_mask, ocean=1, land=0, shape=(89,180)).

xlon,xlat = np.meshgrid(lons,lats)
PP.pcolormesh(xlon,xlat,land_mask)
PP.colorbar()
PP.show()

Land Mask

I would like to loop over all the lats and lons and do a calculation for those lat/lon pairs that are over the ocean, and do nothing ie move on to the next lat/lon pair if over land. Some pseudo code:

# loop over lats
for lat in lats:
    # loop over lons
    for lon in lons:
        # check to see if the current 
        # lat/lon is in the ocean.
        # if True, do something
        if (lat,lon) in ocean:
            do something
        # else if current lat/lon
        # over land, do nothing and
        # move to the next lat/lon
        # pair
        else: # ie if over land, skip this point
            continue

Im not sure how to do this with the 2d land mask that I have. Also perhaps there is a better way to implement this code that is faster and/or more pythonic than nested for loops? Thanks in advance.

2
pythonnumpymultidimensional-arraylogical-operators

2 Answers

1
votes

I imagine an idea like below would work.

import numpy
a = numpy.arange(9).reshape(3, 3)
# array([[0, 1, 2],
#        [3, 4, 5],
#        [6, 7, 8]])
b = a > 4
# array([[False, False, False],
#        [False, False,  True],
#        [ True,  True,  True]])
c = numpy.zeros(a.shape)
# array([[0., 0., 0.],
#        [0., 0., 0.],
#        [0., 0., 0.]])
c[~b] = 1
# array([[1., 1., 1.],
#        [1., 1., 0.],
#        [0., 0., 0.]])

This means that you can use your mask to modify only specific entries of an array which has the same shape.

0
votes

You could try this:

nlats, nlons = land_mask.shape
for i in range(nlons):
   for j in range(nlats):
       if land_mask(j,i) == 1:
           do something

But this will be very slow in python. Numpy operations can be made faster by replacing loops with vectorization. So, there might be a better way to do this if you tell us what the do something part is supposed to do.