Python Cartopy mapping

Question

Я пытаюсь отобразить уровни CO2 на основе общедоступных данных NASA на глобальной карте и изобразить эти значения как (long, lat, value) как топографическую карту, основываясь на данных и с помощью программного обеспечения Panoply, это то, что мой график долженвыглядит так:

Данные в формате .nc4 и читаются правильно, однако я не могу получить график данных Я использую Cartopy API и следую этому примеру: (https://scitools.org.uk/cartopy/docs/latest/gallery/waves.html#sphx-glr-gallery-waves-py).

Также я не хочу использовать Базовую карту.

Попытка # 1:

См. Код Python ниже:

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs


"""
function that download each OCO - 2 data that is in .nc4 format from file "subset_OCO2_L2_ABand_V8_20180929_010345.txt"
which is list of links
for all data with date range 2015 - 09 - 01 to 2016 - 01 - 01# make sure that you have a valid user name & password by registering in https: //earthdata.nasa.gov/
#implementation based on http: //unidata.github.io/netcdf4-python/#section1"""

def download_oco2_nc4(username, password, filespath):
  filespath = "C:\\Users\\Desktop\\oco2\\oco2_LtCO2_150831_B8100r_171009083146s.nc4"
  dataset = Dataset(filespath)
  print(dataset.file_format)
  print(dataset.dimensions.keys())
  print(dataset.variables['xco2'])
  XCO2 = []
  LONGITUDE = []
  LATITUDE = []
  # XCO2
  XCO2 = dataset.variables['xco2'][:]
  print("->", type(XCO2))
  print(dataset.variables['latitude'])
  # LATITUDE
  LATITUDE = dataset.variables['latitude'][:]
  print(dataset.variables['longitude'])
  # LONGITUDE
  LONGITUDE = dataset.variables['longitude'][:]
  return XCO2, LONGITUDE, LATITUDE, dataset


def mapXoco2():
  fig = plt.figure(figsize = (10, 5))
  ax = fig.add_subplot(1, 1, 1, projection = ccrs.Mollweide())

  XCO2, LONGITUDE, LATITUDE, dataset = download_oco2_nc4(1, 2, 3)
  dataset.close()

  XCO2_subset = list()
  counter = 0
  for xco2 in XCO2:
      if counter < 10:
          XCO2_subset.append(xco2)
          counter = counter + 1
      else:
          break
  print("XCO2_subset="+str(len(XCO2_subset)))


  counter = 0
  LONGITUDE_subset = list()
  for longitude in LONGITUDE:
      if counter < 10:
          LONGITUDE_subset.append(longitude)
          counter = counter + 1
      else:
          break
  print("LONGITUDE_subset="+str(len(LONGITUDE_subset)))

  counter = 0
  LATITUDE_subset = list()
  for latitude in LATITUDE:
      if counter < 10:
          LATITUDE_subset.append(latitude)
          counter = counter + 1
      else:
          break
  print("LATITUDE_subset="+str(len(LATITUDE_subset)))

  XCO2_subset = np.array(XCO2_subset)
  LONGITUDE_subset = np.array(LONGITUDE_subset)
  LATITUDE_subset = np.array(LATITUDE_subset)
  #LONGITUDE_subset, LATITUDE_subset = np.meshgrid(LONGITUDE_subset, LATITUDE_subset)
  #XCO2_subset,XCO2_subset = np.meshgrid(XCO2_subset,XCO2_subset)
  ax.contourf(LONGITUDE_subset,LATITUDE_subset,XCO2_subset,
      transform = ccrs.Mollweide(central_longitude=0, globe=None),
      cmap = 'nipy_spectral')
  ax.coastlines()
  ax.set_global()
  plt.show()
  print(XCO2_subset)

mapXoco2()

Когда я закомментирую эти строки:

 #LONGITUDE_subset, LATITUDE_subset = np.meshgrid(LONGITUDE_subset, LATITUDE_subset)
      #XCO2_subset,XCO2_subset = np.meshgrid(XCO2_subset,XCO2_subset)

Я получаю ошибку:

повышение TypeError («Input z должен быть 2D-массивом.»)

TypeError: Input z должно быть 2D-массивом.

Однакокогда я НЕ комментирую эти строки:

 LONGITUDE_subset, LATITUDE_subset = np.meshgrid(LONGITUDE_subset, LATITUDE_subset)
      XCO2_subset,XCO2_subset = np.meshgrid(XCO2_subset,XCO2_subset

)

Я получаю пустую карту, я вижу континенты, но нет отложенных значений Значения C02.

Я считаю, что неправильно интерпретировал 1D в 2D преобразование моих входов.

Попытка № 2 (обновлено):

Вместо того, чтобы понимать, почему / что делают эти 2d-преобразования в API, я строю каждую точку 1 на 1, используя цикл.Проблема в том, что я могу видеть больше данных (я строю только около 10% данных) Я НЕ МОГУ ВИДЕТЬ КАРТУ / КОНТИНЕНТЫ Я ВИДЮ ЗНАЧЕНИЕ ЗНАЧЕНИЯ НА БЕЛОМ ФОНЕ ???* from netCDF4 import Dataset import numpy as np import matplotlib.pyplot as plt import numpy as np import cartopy.crs as ccrs from random import sample """ function that download each OCO - 2 data that is in .nc4 format from file "subset_OCO2_L2_ABand_V8_20180929_010345.txt" which is list of links for all data with date range 2015 - 09 - 01 to 2016 - 01 - 01# make sure that you have a valid user name & password by registering in https: //earthdata.nasa.gov/ #implementation based on http: //unidata.github.io/netcdf4-python/#section1""" filespath = "C:\\Users\\Downloads\\oco2_LtCO2_150830_B7305Br_160712072205s.nc4" def download_oco2_nc4(filespath): dataset = Dataset(filespath) print("file format:"+str(dataset.file_format)) print("dimensions.keys():"+str(dataset.dimensions.keys())) print("variables['xco2']:"+str(dataset.variables['xco2'])) XCO2 = [] LONGITUDE = [] LATITUDE = [] # XCO2 XCO2 = dataset.variables['xco2'][:] print("->", type(XCO2)) print(dataset.variables['latitude']) # LATITUDE LATITUDE = dataset.variables['latitude'][:] print(dataset.variables['longitude']) # LONGITUDE LONGITUDE = dataset.variables['longitude'][:] return XCO2, LONGITUDE, LATITUDE, dataset def mapXoco2(): fig = plt.figure(figsize = (10, 5)) ax = fig.add_subplot(1, 1, 1, projection = ccrs.Mollweide()) XCO2, LONGITUDE, LATITUDE, dataset = download_oco2_nc4(filespath) dataset.close() XCO2_subset = np.array(XCO2) LONGITUDE_subset = np.array(LONGITUDE) LATITUDE_subset = np.array(LATITUDE) """each of the arrays has over 80,000 of data therefore its taking to long to map, after 10,000 rows its to slow, and 10,000 isnt sufficient. Because oco-2 gathers data from trajectory the 1st 10% or whatever precent of the data will not be a good representation of the overal data. We must sample from X number of slices across the data. """ #XCO2 attempt to get ten ranges, we need to check 10 ranges therefore we need if statements not if/else if (len(XCO2_subset)>=10000): first_XCO2_subset=XCO2_subset[0:1000] if (len(XCO2_subset)>=20000): second_XCO2_subset=XCO2_subset[20000:21000] if (len(XCO2_subset)>=30000): third_XCO2_subset=XCO2_subset[30000:31000] if (len(XCO2_subset)>=40000): fourth_XCO2_subset=XCO2_subset[40000:41000] if (len(XCO2_subset)>=50000): fifth_XCO2_subset=XCO2_subset[50000:51000] if (len(XCO2_subset)>=60000): sixth_XCO2_subset=XCO2_subset[60000:61000] if (len(XCO2_subset)>=70000): seventh_XCO2_subset=XCO2_subset[70000:71000] if (len(XCO2_subset)>=80000): eight_XCO2_subset=XCO2_subset[80000:81000] sampled_xco2 = first_XCO2_subset + second_XCO2_subset + third_XCO2_subset + fourth_XCO2_subset + fifth_XCO2_subset + sixth_XCO2_subset + seventh_XCO2_subset + eight_XCO2_subset #LONGITUDE attempt to get ten ranges, we need to check 10 ranges therefore we need if statements not if/else if (len(LONGITUDE_subset)>=10000): first_LONGITUDE_subset=LONGITUDE_subset[0:1000] if (len(LONGITUDE_subset)>=20000): second_LONGITUDE_subset=LONGITUDE_subset[20000:21000] if (len(LONGITUDE_subset)>=30000): third_LONGITUDE_subset=LONGITUDE_subset[30000:31000] if (len(LONGITUDE_subset)>=40000): fourth_LONGITUDE_subset=LONGITUDE_subset[40000:41000] if (len(LONGITUDE_subset)>=50000): fifth_LONGITUDE_subset=LONGITUDE_subset[50000:51000] if (len(LONGITUDE_subset)>=60000): sixth_LONGITUDE_subset=LONGITUDE_subset[60000:61000] if (len(LONGITUDE_subset)>=70000): seventh_LONGITUDE_subset=LONGITUDE_subset[70000:71000] if (len(LONGITUDE_subset)>=80000): eight_LONGITUDE_subset=LONGITUDE_subset[80000:81000] sampled_LONGITUDE = first_LONGITUDE_subset + second_LONGITUDE_subset + third_LONGITUDE_subset + fourth_LONGITUDE_subset + fifth_LONGITUDE_subset + sixth_LONGITUDE_subset + seventh_LONGITUDE_subset + eight_LONGITUDE_subset #LATITUDE attempt to get ten ranges, we need to check 10 ranges therefore we need if statements not if/else if (len(LATITUDE_subset)>=10000): first_LATITUDE_subset=LATITUDE_subset[0:1000] if (len(LATITUDE_subset)>=20000): second_LATITUDE_subset=LATITUDE_subset[20000:21000] if (len(LATITUDE_subset)>=30000): third_LATITUDE_subset=LATITUDE_subset[30000:31000] if (len(LATITUDE_subset)>=40000): fourth_LATITUDE_subset=LATITUDE_subset[40000:41000] if (len(LATITUDE_subset)>=50000): fifth_LATITUDE_subset=LATITUDE_subset[50000:51000] if (len(LATITUDE_subset)>=60000): sixth_LATITUDE_subset=LATITUDE_subset[60000:61000] if (len(LATITUDE_subset)>=70000): seventh_LATITUDE_subset=LATITUDE_subset[70000:71000] if (len(LATITUDE_subset)>=80000): eight_LATITUDE_subset=LATITUDE_subset[80000:81000] sampled_LATITUDE = first_LATITUDE_subset + second_LATITUDE_subset + third_LATITUDE_subset + fourth_LATITUDE_subset + fifth_LATITUDE_subset + sixth_LATITUDE_subset + seventh_LATITUDE_subset + eight_LATITUDE_subset ax = plt.axes(projection=ccrs.Mollweide()) #plt.contourf(LONGITUDE_subset, LATITUDE_subset, XCO2_subset, 60,transform=ccrs.PlateCarree()) for long, lat, value in zip(sampled_LONGITUDE, sampled_LATITUDE,sampled_xco2): #print(long, lat, value) if value >= 0 and value < 370: ax.plot(long,lat,marker='o',color='blue', markersize=1, transform=ccrs.PlateCarree()) elif value >= 370 and value < 390: ax.plot(long,lat,marker='o',color='cyan', markersize=1, transform=ccrs.PlateCarree()) elif value >= 390 and value < 402: ax.plot(long,lat,marker='o',color='yellow', markersize=1, transform=ccrs.PlateCarree()) elif value >= 402 and value < 410: ax.plot(long,lat,marker='o',color='orange', markersize=1, transform=ccrs.PlateCarree()) elif value >= 410 and value < 415: ax.plot(long,lat,marker='o',color='red', markersize=1, transform=ccrs.PlateCarree()) else: ax.plot(long,lat,marker='o',color='brown', markersize=1, transform=ccrs.PlateCarree()) ax.coastlines() plt.show() mapXoco2() Вывод:

формат файла: NETCDF4

sizes.keys (): odict_keys (['sounding_id', 'levels', 'band',' vertices ',' epoch_dimension ',' source_files '])

переменные [' xco2 ']: float32 xco2 (sounding_id) единицы: ppm long_name: XCO2 missing_value: -999999.0 комментарий: усредненный по столбцам dry-воздушная мольная доля CO2 (включая коррекцию смещения)

неограниченные размеры: текущая форма = (82776,) заполнение включено, значение по умолчанию _FillValue 9,969209968386869e + 36 используется

-> широта float32 (sounding_id)Единицы измерения: градусы_северное длинное_имя: широта отсутствует_значение: -999999.0 комментарий: центральная широта измерения неограниченные размеры: текущая форма = (82776,) заполнение включено, значение по умолчанию _FillValue 9,969209968386869e + 36 используемых

float32 долгота (sounding_id) единиц: degrees_east long_name: longitude missing_value: -999999.0 комментарий: центральная долгота измерения неограниченные размеры: текущая форма = (82776,) заполнение включено, по умолчанию _FillValue 9,969209968386869e + 36 используется

1) Что случилось с картой и континентами?

Спасибо и любая полезная помощь оценена.

Answer 1

Похоже, ваш аргумент преобразования неверен.Если у вас есть данные широты / долготы, значение аргумента преобразования должно быть ccrs.PlateCarree().См. Это руководство в документации по картоплате для деталей: https://scitools.org.uk/cartopy/docs/latest/tutorials/understanding_transform.html.

Я не могу запустить ваш пример, чтобы убедиться, что это решение.Чтобы получить максимальную отдачу от переполнения стека, вы должны предоставить минимальный рабочий пример, чтобы другие могли работать сами.См. https://stackoverflow.com/help/mcve и https://stackoverflow.com/help/how-to-ask для подсказок.