Im trying to project Goes-16 Netcdf files into basemap, but im running into the same error every time. Ive tried downgrading python from 3.7.3 to 3.6, updating and downgrading basemap, and there's no literature of this matter anywhere I could find.
I read online that maybe python 3.7.3 might not be compatible with basemap anymore, but I ran the code with python 3.7.1 before and it worked. I tried looking for the alternative mentioned in the code, which is "inspect.cleandoc", but I can't figure out how to properly use this library.
My code is:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
from remap import remap
from cpt_convert import loadCPT
from matplotlib.colors import LinearSegmentedColormap
from netCDF4 import Dataset
from matplotlib.patches import Rectangle
from osgeo import gdal
# Load the Data
# Path to the GOES-16 image file
channel_8 = 'C08/OR_ABI-L2-CMIPF-M3C08_G16_s20180391900384_e20180391911151_c20180391911220.nc'
path_30 = '/Volumes/Anthonys_backup/Masters_Thesis/Satellite_Data/new/'
path = path_30 + channel_8
High_q = 'yes'
geotif = 'no'
"Check remap.py!!!!!!!!"
nc = Dataset(path, 'r',)
H = nc.variables['goes_imager_projection'].perspective_point_height
#band = nc.variables['band_id'][:][0]
#Band = np.char.mod('%d', band)
Band = 8
# Choose the visualization extent (min lon, min lat, max lon, max lat)
extent = [-76, -60, -51, -20]
# Choose the image resolution (the higher the number the faster the processing is)
resolution = .1
# Call the reprojection funcion
grid = remap(path, extent, resolution, 'HDF5')
# Read the data returned by the function
if Band >= 7:
data = grid.ReadAsArray()
Unit = "Brightness Temperature [K]"
else:
data = grid.ReadAsArray()
a = data
# set any negative values to zero\n"
data[data < 0.] = 0.
# normalize data to 1.16
data = data/1.16
# set any values greater than 1.0 equal to 1.0
data[a>1.]=1.
Unit = "Reflectance"
# Define the size of the saved picture
if High_q == 'yes':
DPI = 150
fig = plt.figure(figsize=(data.shape[1]/float(DPI), data.shape[0]/float(DPI)),
frameon=False, dpi=DPI)
ax = plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_off()
fig.add_axes(ax)
ax = plt.axis('off')
else:
DPI = 150
# Plot the Data
# Create the basemap reference for the Rectangular Projection
bmap = Basemap(resolution='h', llcrnrlon=extent[0], llcrnrlat=extent[1], urcrnrlon=extent[2], urcrnrlat=extent[3], epsg=4326)
# Draw the countries and Argentinian states shapefiles
bmap.readshapefile('/Users/anthonycrespo/Desktop/New_code_for_plotting/arg_adm1/ARG_adm1','ARG_adm1',linewidth=.5,color='black')
bmap.readshapefile('/Users/anthonycrespo/Desktop/New_code_for_plotting/Countries_Shape/ne_10m_admin_0_countries','ne_10m_admin_0_countries',linewidth=.7,color='black')
# Draw parallels and meridians
bmap.drawcoastlines(linewidth=1., linestyle='solid', color='black')
bmap.drawparallels(np.arange(-90.0, 90.0, 3.), linewidth=0.3,
dashes=[4, 4], color='white', labels=[True,
True, False, False], fmt='%g', labelstyle="+/-",
xoffset=0.10, yoffset=-1.00, size=5)
bmap.drawmeridians(np.arange(0.0, 360.0, 3.), linewidth=0.3,
dashes=[4, 4], color='white', labels=[False,
False, True, False], fmt='%g', labelstyle="+/-",
xoffset=-0.80, yoffset=0.20, size=5)
# Converts a CPT file to be used in Python
if 7<= Band <= 10:
cpt = loadCPT('/Users/anthonycrespo/Desktop/New_code_for_plotting/WVCOLOR35.cpt')
elif Band >= 11:
cpt = loadCPT('/Users/anthonycrespo/Desktop/New_code_for_plotting/IR4AVHRR6.cpt')
else:
pass
# Makes a linear interpolation
if Band >= 7:
cpt_convert = LinearSegmentedColormap('cpt', cpt)
else:
pass
# Plot the GOES-16 channel with the converted CPT colors
if Band >= 8:
bmap.imshow(data, origin='upper', cmap=cpt_convert, vmin=170, vmax=378)
else:
bmap.imshow(data, origin='upper', cmap='Greys_r', vmin=0., vmax=1.)
# Date and time
import datetime
time_var = nc.time_coverage_start
iyear = time_var[0:4]
imonth = time_var[5:7]
import calendar
cmonth = calendar.month_abbr[int(imonth)]
iday = time_var[8:10]
itime = time_var[11:19]
itimehr = time_var[11:13]
itimemn = time_var[14:16]
ctime_string = iyear +' '+cmonth+' '+iday+' '+itime+' UTC'
ctime_file_string = iyear + imonth + iday + itimehr + itimemn
filestring = 'Channel 8' + iyear + imonth + iday + "_" + itimehr + itimemn + ".png"
filestring2 = 'CHANNEL 8' + iyear + imonth + iday + "_" + itimehr + itimemn + ".tif"
time_string = 'GOES-16 ABI Band 8 %s '%ctime_string
name = "GOES-16 ABI Band 8"
Name = (name + ' ' + ' ' + ' ' + ' ' + Unit)
time_stamp = ("%s" %ctime_string)
# Add a black rectangle in the bottom to insert the image description
lon_difference = (extent[2] - extent[0]) # Max Lon - Min Lon
currentAxis = plt.gca()
currentAxis.add_patch(Rectangle((extent[0], extent[1]), lon_difference,
lon_difference * 0.060, alpha=1, zorder=3,
facecolor='black'))
# Add the image description inside the black rectangle
lat_difference = (extent[3] - extent[1]) # Max lat - Min lat
plt.text(extent[0], extent[1] + lat_difference * 0.005, Name,
horizontalalignment='left', color = 'white', size=3)
plt.text(extent[2], extent[1] + lat_difference * 0.005, time_stamp,
horizontalalignment='right', color = 'white', size=3)
# Insert the colorbar at the right
cb = bmap.colorbar(location='bottom', size = '2%', pad = '-4%')
cb.outline.set_visible(False) # Remove the colorbar outline
cb.ax.tick_params(width = 0) # Remove the colorbar ticks
cb.ax.xaxis.set_tick_params(pad=-7.5) # Put the colobar labels inside the colorbar
cb.ax.tick_params(axis='x', colors='grey', labelsize=3) # Change the color and size of the colorbar labels
# Show the plot
#plt.show()
plt.savefig('/Users/anthonycrespo/Desktop/New_code_for_plotting/Images/' + filestring, dpi=DPI, bbox_inches='tight', pad_inches=0)
# Export the result to GeoTIFF
if geotif == 'yes':
driver = gdal.GetDriverByName('GTiff')
driver.CreateCopy('/Users/anthonycrespo/Desktop/New_code_for_plotting/Geotif/' + filestring2, grid, 0)
else:
pass
# Close the plot window
plt.close()
And the complete error message is:
ERROR 5: GDALWarpOptions.Validate(): nBandCount=0, no bands configured!
- finished! Time: 0.00048089027404785156 seconds
Reprojection.py:75: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
bmap = Basemap(resolution='h', llcrnrlon=extent[0], llcrnrlat=extent[1], urcrnrlon=extent[2], urcrnrlat=extent[3], epsg=4326)
Reprojection.py:78: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
bmap.readshapefile('/Users/anthonycrespo/Desktop/New_code_for_plotting/arg_adm1/ARG_adm1','ARG_adm1',linewidth=.5,color='black')
Reprojection.py:79: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
bmap.readshapefile('/Users/anthonycrespo/Desktop/New_code_for_plotting/Countries_Shape/ne_10m_admin_0_countries','ne_10m_admin_0_countries',linewidth=.7,color='black')
Segmentation fault: 11