1. I would like to know how I can create a gridded map of a country(i.e. Singapore) with resolution of 200m x 200m squares. (50m or 100m is ok too)

2. I would then use the 'nearest neighbour' technique to assign a rainfall data and colour code to each square based on the nearest rainfall station's data. [I have the latitude,longitude & rainfall data for all the stations for each date.]

3. Then, I would like to store the data in an Array for each 'gridded map' (i.e. from 1-Jan-1980 to 31-Dec-2021)

Can this be done using python?

P.S Below is a 'simple' version I did as an example to how the 'gridded' map should look like for 1 particular day.

data sourcing

import requests, itertools, io
from pathlib import Path
import urllib
from zipfile import ZipFile
import fiona.drvsupport
import geopandas as gpd
import numpy as np
import pandas as pd
import shapely.geometry

# get official Singapore planning area geometry
url = "https://geo.data.gov.sg/planning-area-census2010/2014/04/14/kml/planning-area-census2010.zip"

f = Path.cwd().joinpath(urllib.parse.urlparse(url).path.split("/")[-1])
if not f.exists():
    r = requests.get(url, stream=True, headers={"User-Agent": "XY"})
    with open(f, "wb") as fd:
        for chunk in r.iter_content(chunk_size=128):
            fd.write(chunk)
zfile = ZipFile(f)
zfile.extractall(f.stem)

fiona.drvsupport.supported_drivers['KML'] = 'rw'
gdf_sg = gpd.read_file(
    [_ for _ in Path.cwd().joinpath(f.stem).glob("*.kml")][0], driver="KML"
)


# get data about Singapore weather stations
df_stations = pd.json_normalize(
    requests.get("https://api.data.gov.sg/v1/environment/rainfall").json()["metadata"][
        "stations"
    ]
)

# dates to get data from weather.gov.sg
dates = pd.date_range("20220601", "20220730", freq="MS").strftime("%Y%m")
df = pd.DataFrame()
# fmt: off
bad = ['S100', 'S201', 'S202', 'S203', 'S204', 'S205', 'S207', 'S208',
       'S209', 'S211', 'S212', 'S213', 'S214', 'S215', 'S216', 'S217',
       'S218', 'S219', 'S220', 'S221', 'S222', 'S223', 'S224', 'S226',
       'S227', 'S228', 'S229', 'S230', 'S900']
# fmt: on
for stat, month in itertools.product(df_stations["id"], dates):
    if not stat in bad:
        try:
            df_ = pd.read_csv(
                io.StringIO(
                    requests.get(
                        f"http://www.weather.gov.sg/files/dailydata/DAILYDATA_{stat}_{month}.csv"
                    ).text
                )
            ).iloc[:, 0:5]
        except pd.errors.ParserError as e:
            bad.append(stat)
            print(f"failed {stat} {month}")
        df = pd.concat([df, df_.assign(id=stat)])

df["Rainfall (mm)"] = pd.to_numeric(
    df["Daily Rainfall Total (mm)"], errors="coerce"
)
df["Date"] = pd.to_datetime(df[["Year","Month","Day"]]).dt.strftime("%Y%m%d")
df = df.loc[:,["id","Date","Rainfall (mm)", "Station"]]

visualisation using plotly animation

import plotly.express as px

# reduce dates so figure builds in sensible time
gdf_px = gdf_grid_rainfall.loc[
    lambda d: d["Date"].isin(
        gdf_grid_rainfall["Date"].value_counts().sort_index().index[0:15]
    )
]

px.choropleth_mapbox(
    gdf_px,
    geojson=gdf_px.geometry,
    locations=gdf_px.index,
    color="Rainfall (mm)",
    hover_data=gdf_px.columns[1:].tolist(),
    animation_frame="Date",
    mapbox_style="carto-positron",
    center={"lat":gdf_px.unary_union.centroid.y, "lon":gdf_px.unary_union.centroid.x},
    zoom=8.5
).update_layout(margin={"r": 0, "t": 0, "l": 0, "b": 0, "pad": 4})