1

I have an issue when adapting a simple dataset to generate temperature interpolation map from points


#!apt-get update
#!apt-get -qq install python-cartopy python3-cartopy

#!pip install metpy

installed dependecies as using a colab notebook

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

import requests
from pandas.io.json import json_normalize
import pandas as pd

!curl "https://climaya.com/fp-weather/api/api.php?q=GT_forecast&d=3" > data.json
df = pd.read_json("data.json")
#df

lat=df['latitude'].values
lon=df['longitude'].values
#lat,lon

this gives an array of lat, lon

mapcrs=ccrs.PlateCarree(central_longitude=-90)
datacrs = ccrs.PlateCarree()

fig=plt.figure(figsize=(10,8))
ax=fig.add_subplot(1,1,1,projection=mapcrs)
ax.scatter(lon,lat,c=df['temp_high'].values*1,transform=datacrs)

this shows the data points ok

enter image description here


from metpy.interpolate import interpolate_to_grid, remove_nan_observations

lon=df['longitude'].values*1.
lat=df['latitude'].values*1.
dat=df['temp_high'].values*1.0

#lon,lat,dat data seems ok shape is 56, each just to asure numeric values

xp, yp, _ = mapcrs.transform_points(datacrs, lon, lat).T

x_masked, y_masked, dat = remove_nan_observations(xp, yp, dat)

gridx, gridy, gridz = interpolate_to_grid(x_masked, y_masked, dat, interp_type='cressman', minimum_neighbors=1, search_radius=40000, hres=5000)
#command in notebook ok no errors

import numpy as np
fig=plt.figure(figsize=(10,8))
ax=fig.add_subplot(1,1,1,projection=mapcrs)
ax.contourf(gridx,gridy,gridz,np.arange(0,35.1,1))

this generate an error TypeError: Input z must be at least a (2, 2) shaped array, but has shape (1, 1)

as new to python and metpy have struggled with the solution, thanks for any help.

**EDITED: Thanks @dopplerShift, that was the thing, here is the fix and the result

gridx, gridy, gridz = interpolate_to_grid(lon, lat, dat, interp_type='cressman', minimum_neighbors=1, search_radius=.75, hres=.1)

enter image description here

neavilag
  • 609
  • 1
  • 8
  • 20

1 Answers1

1

So the problem is this line from the docs:

hres (float) – The horizontal resolution of the generated grid, given in the same units as the x and y parameters. Default 50000.

Here you pass hres=5000, which is saying "generate a grid with 5000 units between points". For points coming through a map projection, this is generally meters, except that you're using PlateCarree(), where the points are given in degrees. So for your current map projection, you should use an hres value in degrees, so maybe 0.25 or 0.1. You'll also want to adjust search_radius to a value in degrees (rather than the current 40km) as well--or leave it out and it will try to calculate something sensible.

DopplerShift
  • 5,472
  • 1
  • 21
  • 20