Create a table of longitudes and latitudes and implement the cells corresponding to the coordinates within a country (python)

Question

I could use some help. I have created a two-dimensional array filled with 0 which has as dimension the longitudes and latitudes of the northern hemisphere with a resolution of 0.5°. My goal is to go through the latitudes and longitudes of this array so that when they are included in the boundaries of France, the 0 is changed to 1. The final goal is to create a NetCdf file with the longitudes, latitudes and binary values as columns. The data of this file would then be multiplied by an inventory to have this inventory only for France.

The code I wrote for now is the one below but it seems to have a problem with the polygon because when I try with only 4 sets of coordinates, it works.

Thanks in advance for your help

import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import Polygon, Point, MultiPolygon
import shapefile

lon = np.arange(-179.75,180.25,0.5)
lat = np.arange(0.25,89.25,0.5)
tab = np.zeros((len(lon), len(lat))) #creation of the array filled with 0 with the dimensions len(lon) and len(lat)
 
shapefile = gpd.read_file("./TM_WORLD_BORDERS-0.3.shp") #open the shapefile containing the borders of the countries
#print(shapefile)

def SelectCountry(Country) : #function to create a dataset with only the desired country 
    mask = (shapefile['NAME'] == Country)
    p = shapefile[mask]
    return p

france = SelectCountry('France') #selection of France 

for i in range(len(lon)) : #Put 1 in the cells of the table corresponding to the latitudes and longitudes included in the shape of France
    for j in range(len(lat)):
        point = shapely.geometry.Point(lon[i], lat[j])
        shape = france['geometry']
        if shape.contains(point):
            tab[i,j] = 1
                
np.where(tab == 1)
print(tab)

Answer 1

check out regionmask - it's designed to integrate with xarray (which you should also take a look at) to facilitate transitioning from gridded to polygon-based data definitions.

If you want to just use numpy arrays, you could use shapely.vectorized.contains to build your boolean mask:

XX, YY = np.meshgrid(lon, lat)
in_polygon = shapely.vectorized.contains(france['geometry'].item(), XX, YY)

in_polygon # boolean array with shape (Nlons, Nlats)