Почему гексбиновая сетка искажается в полярной стереопроекции - PullRequest
Новая
       41

Почему гексбиновая сетка искажается в полярной стереопроекции

0 голосов
/ 09 января 2019

Я пытаюсь понять, почему на гексбинном графике в северной или южной полярной стереопроекции отображаются сдавленные шестиугольники, хотя площадь сетки квадратная, а проекция примерно равна площади.

Я попробовал стереопроекции как в северной, так и в южной полярности, используя базовую карту. 

import numpy as np
from numpy.random import uniform
import matplotlib.pyplot as plt 
from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(12,10)) # width, height in inches
ax =fig.add_axes([-0.02,0.1,0.74,0.74]) 

m = Basemap(epsg='3413',lon_0=0.,resolution='l',width=6000000,height=6000000)

m.drawcoastlines()

m.drawmapscale(0.,90.,0.,90.,1000)

npts=2000
lats = uniform(60.,80.,size=npts)
lons = uniform(0.,360.,size=npts)
data = uniform(0.,4800.,size=npts)

x,y=m(lons, lats)

thiscmap=plt.cm.get_cmap('viridis')

p=m.hexbin(x,y,C=data,gridsize=[10,10],cmap=thiscmap)

plt.show()

plot output

1 Ответ

0 голосов
/ 10 января 2019

Я не знаю, почему вы получаете squashed hexagons. Но вы можете настроить форму шестиугольника, установив соответствующие значения gridsize. Здесь я изменяю ваш код и получаю лучший сюжет. 

import numpy as np
from numpy.random import uniform
import matplotlib.pyplot as plt 
from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(12,10)) # width, height in inches
ax =fig.add_axes([-0.02, 0.1, 0.74, 0.74]) 

# North polar stereographic projection epsg='3413'; ***large areal distortion***
#m = Basemap(epsg='3413', lon_0=0., resolution='c', width=6000000, height=6000000)

# 'laea':  Lambert Azimuthal Equal Area
# Thematic mapping with ground surface data should be plotted on 'equal-area' projection
m = Basemap(projection='laea', lon_0=0., lat_0=90, resolution='l', width=6000000, height=6000000)

m.drawcoastlines(linewidth=0.5)

m.drawmapscale(0.,90.,0.,90.,1000)  # 1000 km?

npts = 2000
lats = uniform(60.,80.,size=npts)  # not cover N pole
lons = uniform(0.,360.,size=npts)  # around W to E
data = uniform(0.,4800.,size=npts)

x,y = m(lons, lats)

thiscmap = plt.cm.get_cmap('viridis')

# To get 'rounded' hexagons, gridsize should be specified appropriately
# need some trial and error to get them right
#p=m.hexbin(x, y, C=data, gridsize=[10,10], cmap=thiscmap)  # original code
m.hexbin(x, y, C=data, gridsize=[16,11], cmap=thiscmap)     # better

plt.colorbar()  # useful on thematic map

plt.show()

Используемая вами проекция (epsg: 3413) - это стереографическая проекция с большими ареалами ареала. Более подходящая проекция для тематического картирования в этом случае - Lambert Azimuthal Equal Area.

Полученный участок:

enter image description here

...