Thanks! That did it. The problem was:
# specify model
with pm.Model(coords=coords) as sk_lm:
# data
feature_data = pm.Data('feature_data', features, dims=('obs_id', 'features'))
# priors
alpha = pm.Normal('alpha', mu=0, sigma=1)
betas = pm.Normal('betas', mu=[40, 90, 50], sigma=5, dims='features')
# model error
sigma = pm.Exponential("sigma", lam=1)
# matrix-dot products
m1 = pm.math.matrix_dot(feature_data, betas) # here
I was using the literal features data instead of the pm.Data object.