1

Exactly as the title says. I'm plotting sea ice concentration data and this code:

map = Basemap(projection='cyl', lat_0 = lat_0, lon_0 = lon_0,
          llcrnrlon = llcrnrlon, llcrnrlat = llcrnrlat,
          urcrnrlat = urcrnrlat, urcrnrlon = urcrnrlon,
          area_thresh = 1000., resolution='l')

Works just fine. https://i.stack.imgur.com/9RcSC.png

But when I try to change the projection:

map = Basemap(projection='npstere',boundinglat=10,lon_0=270,resolution='l')

I just get a blank map. https://i.stack.imgur.com/s42Q7.png

Am I forgetting something about basemap? I though changing the projection was pretty simple.

Edit: Heres the full code

import numpy as np
import math as m
import urllib2
import time
import datetime
import matplotlib.pyplot as plt
import matplotlib.patheffects as path_effects
from matplotlib.colors import LinearSegmentedColormap
from mpl_toolkits.basemap import Basemap
from pydap.client import open_url
from pydap.proxy import ArrayProxy
import scipy
from scipy.ndimage.filters import minimum_filter, maximum_filter

data_url_ice = 'http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/noaa.oisst.v2.highres/icec.day.mean.2015.v2.nc'

dataset3 = open_url(data_url_ice)

#############################################################################################
# Get Data
#############################################################################################

# Daily Mean Ice Concentration

Lat_ice = dataset3['lat']
Lon_ice = dataset3['lon']
Time_ice = dataset3['time']
Ice_Conc = dataset3['icec']
Ice_Conc = Ice_Conc[-1,:,:]
Ice_Conc = Ice_Conc.array[:]

Ice_Conc[Ice_Conc < 0] = 0
Ice_Conc = Ice_Conc * 100.

Ice_Conc = Ice_Conc.squeeze()

#############################################################################################
# Colormap
#############################################################################################

vmax_ice = 100.0
cmap_ice = LinearSegmentedColormap.from_list('mycmap', [(0     / vmax_ice, 'white'),    #-40
                                                         (50.   / vmax_ice, 'yellow'), #-20
                                                         (100.   / vmax_ice, 'blue')]    # 20
                                              )

#############################################################################################
# Map Projection Info
#############################################################################################

lat_0 = 0
lon_0 = 0

llcrnrlat = -90. # (1,1)
llcrnrlon = 0. # (1,1)

urcrnrlat = 90. # (720,361)
urcrnrlon = 359. # (720,361)

# Daily Mean Ice Concentration

fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1,1,1)
#map = Basemap(projection='ortho', lat_0 = 50, lon_0 = -105,
#              area_thresh = 1000., resolution='i')

map = Basemap(projection='npstere',boundinglat=10,lon_0=270,resolution='l')

#map = Basemap(projection='cyl', lat_0 = lat_0, lon_0 = lon_0,
#              llcrnrlon = llcrnrlon, llcrnrlat = llcrnrlat,
#              urcrnrlat = urcrnrlat, urcrnrlon = urcrnrlon,
#              area_thresh = 1000., resolution='l')

map.drawcoastlines()
map.drawcountries()
map.drawmapboundary()
map.fillcontinents(color='white')

levels = np.linspace(0,100,100)
ticks = [0,10,20,30,40,50,60,70,80,90,100]
iceconc = plt.contourf(Lon_ice,Lat_ice,Ice_Conc,levels,cmap=cmap_ice)

# Set Colorbar Text Color
color_bar = map.colorbar(iceconc)
color_bar.set_ticks(ticks)
cbytick_obj = plt.getp(color_bar.ax.axes, 'yticklabels')
plt.setp(cbytick_obj, color='w')

ax.text(0,1.02,'Sea Ice Concentration (%)\n\n',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax.transAxes, color='w', fontsize=9).set_clip_on(False)
ax.text(0,1.02,'Satellite Measured Daily Mean\n',
        verticalalignment='bottom', horizontalalignment='left',
        transform=ax.transAxes, color='w', fontsize=9).set_clip_on(False)        

plt.savefig('/home/ice_current.png', facecolor='#3e3e3e', bbox_inches='tight')
plt.close("all")
pythonismyjam
  • 99
  • 1
  • 10

1 Answers1

2

I think you're forgetting the simplest line, however you haven't shown most of your code so it's hard to guess. i.e.

map = Basemap(projection='cyl', lat_0 = lat_0, lon_0 = lon_0,
          llcrnrlon = llcrnrlon, llcrnrlat = llcrnrlat,
          urcrnrlat = urcrnrlat, urcrnrlon = urcrnrlon,
          area_thresh = 1000., resolution='l')

#valid aproximate coordinates for some map projections
lat = [19.2325]
lon = [-155.3395]

#but the actual map coordinates are
x,y = map(lat, lon)
map.plot(x, y, 'ro', markersize=6)

I believe your data is still there, but hasn't been converted to proper map coordinates. Notice how in the example they use m, your instance of Basemap is named map in your example so use that. I think it's pretty cool how they did that.

ljetibo
  • 3,048
  • 19
  • 25