Related to my comment from before, maybe we should also return ``NaN`` outside of bounds and make sure that the fitting works with ``NaN``. I'm not saying that this should be done as part of this PR, since it may require changes in ``astropy.modeling``. @adonath - what do you think? Could we somehow adapt ``fit`` below to make it so that NaN values in the data and model are ignored? We could fit just the pixels where both the data and model are not NaN?