Просто могло бы быть, если бы территория была локальной, одно море например: тогда с приемлемой точностью можно было одной командой перепроецировать в локальную проекцию, и брать декартову площадь в ней, умножив количество пикселей попавших в маску на площадь одного пикселя.
Из того, что будет работать, но придётся потратить некоторое время и сделать себе код:
1. numpy 2d массив -> GDAL датасет (не проблема, т.к. пространственные параметры матрицы известны)
2. gdal_polygonize на результат п.1, получается векторный слой, где объекты под маской имеют значение атрибута, равное значениям пикселей в этой маске.
3. через python-ogr открыть вектор из п.2, удалить все объекты которые по значению атрибута не попали в маску
4. В цикле пробегать для всех исходных nympy-матриц пункты 1-3, итоговые векторные объекты закидывая в новый (единый) векторный слой, помечая в атрибутах принадлежность объекта матрице (дата, или идентфикатор какой).
5. Открыть получившийся в п.4 векторный слой в QGIS, настроить вычисление площадей на эллипсоиде, рассчитать для всех объектов площадь через калькулятор полей ($area).
Может кто подскажет, умеет ли pyproj или другая библиотека в сфероидическую геодезию над объектами-полигонами, или геометриями доставаемыми OGR-ом. Если есть что-то - то 4,5 не нужны, сразу вектор из 3 в соответствующую функцию и вот она площадь.