(base) zeng@zeng-virtual-machine:~/soft/pyroms_soft/pyroms/MERRA-2$ python3 get_MERRA_albedo_from_nasa_opendap_daily.py 2022
url https://goldsmr4.gesdisc.eosdis.nasa.go ... 220101.nc4
Traceback (most recent call last):
File "/home/zeng/soft/pyroms_soft/pyroms/MERRA-2/get_MERRA_albedo_from_nasa_opendap_daily.py", line 57, in <module>
lon = dataset['lon'][:]
File "/home/zeng/anaconda3/lib/python3.9/site-packages/pydap/model.py", line 320, in __getitem__
out.data = self._get_data_index(index)
File "/home/zeng/anaconda3/lib/python3.9/site-packages/pydap/model.py", line 349, in _get_data_index
return self._data[index]
File "/home/zeng/anaconda3/lib/python3.9/site-packages/pydap/handlers/dap.py", line 142, in __getitem__
raise_for_status(r)
File "/home/zeng/anaconda3/lib/python3.9/site-packages/pydap/net.py", line 33, in raise_for_status
raise HTTPError(
webob.exc.HTTPError: 302 Found
<!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN">
<html><head>
<title>302 Found</title>
</head><body>
<h1>Found</h1>
<p>The document has moved <a href="https://urs.earthdata.nasa.gov/oauth/au ... re</a>.</p>
</body></html>
I guess there is something wrong with URL?? After I logged into the website, I found that the DATA URL (https://goldsmr4.gesdisc.eosdis.nasa.go ... 220101.nc4) is consistent with the URL created by the code.. So I cannot find the problems, is there anyone willing to help me? Thank you sooooo much!
The code is shown below:
Code: Select all
import matplotlib
matplotlib.use('Agg')
import numpy as np
import netCDF4
from datetime import datetime
import pyroms
import pyroms_toolbox
#import http.cookiejar
#import netrc
#import urllib.request, urllib.error, urllib.parse
#import re
#import pydap.lib
#from pydap.exceptions import ClientError
#import logging
import sys
from pydap.client import open_url
from pydap.cas.urs import setup_session
year = int(sys.argv[1])
invarname = 'ALBEDO'
outvarname = 'albedo'
outtimename = 'albedo_time'
server = 'https://goldsmr4.gesdisc.eosdis.nasa.gov'
leap = year%4
if leap == 0:
daysinmonth = ([31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
daysinyear = 366
else:
daysinmonth = ([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
daysinyear = 365
#if (year <= 1992):
# file_tag = 'MERRA101'
#elif ((year >= 1993) & (year <= 2000)):
# file_tag = 'MERRA201'
#elif ((year >= 2001) & (year <= 2009)):
# file_tag = 'MERRA301'
#elif (year >= 2010):
# file_tag = 'MERRA300'
file_tag = 'MERRA2_400'
#read grid and variable attributes from the first file
year_tag = '%04d' %year
month_tag = '01'
day_tag = '01'
date_tag = year_tag + month_tag + day_tag
url = server + '/opendap/MERRA2/M2T1NXRAD.5.12.4/' + year_tag + '/' + month_tag + '/' + \
file_tag + '.tavg1_2d_rad_Nx.' + date_tag + '.nc4'
print('url',url)
session = setup_session("username","password",url)
dataset = open_url(url, session =session)
lon = dataset['lon'][:]
#shift data between 0 and 360 deg.
gidx = np.where(np.abs(lon) < 1.0e-10)[0][0]
lon = np.asarray(lon)
lon = lon + 180
lat = dataset['lat'][:]
lat = np.asarray(lat)
spval = dataset[invarname].missing_value
units = dataset[invarname].units
long_name = dataset[invarname].long_name
#get data from NASA opendap
for month in range(1-1,1):
nday = 0
#create ROMS forcing file
month_tag = '%02d' %(month+1)
outfile = 'Forcings/MERRA_' + outvarname + '_3hours_' + year_tag + '_' + month_tag + '.nc'
nc = netCDF4.Dataset(outfile, 'w', format='NETCDF3_64BIT')
nc.Author = sys._getframe().f_code.co_name
nc.Created = datetime.now().strftime("%Y-%m-%d %H:%M:%S")
nc.title = 'MERRA-2 dataset. Modern Era Retrospective-analysis'
nc.createDimension('lon', np.size(lon))
nc.createDimension('lat', np.size(lat))
nc.createDimension(outtimename, None)
nc.createVariable('lon', 'f8', ('lon'))
nc.variables['lon'].long_name = 'longitude'
nc.variables['lon'].units = 'degrees_east'
nc.variables['lon'][:] = lon
nc.createVariable('lat', 'f8', ('lat'))
nc.variables['lat'].long_name = 'latitude'
nc.variables['lat'].units = 'degrees_north'
nc.variables['lat'][:] = lat
nc.createVariable(outtimename, 'f8', (outtimename))
nc.variables[outtimename].units = 'days since 1900-01-01 00:00:00'
nc.variables[outtimename].calendar = 'LEAP'
dstart = pyroms_toolbox.date2jday(datetime(year, month+1, 1, 0, 0))
roms_time = np.arange(dstart, dstart+daysinmonth[month], 3./24)
nc.variables[outtimename][:] = roms_time
nc.createVariable(outvarname, 'f', (outtimename, 'lat', 'lon'), fill_value=spval)
nc.variables[outvarname].long_name = long_name
nc.variables[outvarname].units = units
nc.variables[outvarname].coordinates = 'lon lat'
for day in range(1):
# if year == 2010:
# if ((month+1 >= 6) & (month+1 <= 8)):
# file_tag = 'MERRA301'
# else:
# file_tag = 'MERRA300'
day_tag = '%02d' %(day+1)
date_tag = year_tag + month_tag + day_tag
url = server + '/opendap/MERRA2/M2T1NXRAD.5.12.4/' + year_tag + '/' + month_tag + '/' + \
file_tag + '.tavg1_2d_rad_Nx.' + date_tag + '.nc4'
session = setup_session("username","password",url)
dataset = open_url(url, session =session)
var = dataset[invarname][:]
var = np.asarray(var)
#shift data between 0 and 360 deg.
svar = np.zeros(var.shape)
svar[:,:,:len(lon)-gidx] = var[:,:,gidx:]
svar[:,:,len(lon)-gidx:] = var[:,:,:gidx]
mask = np.invert(np.ma.masked_values(svar, spval).mask).astype('float') #mask
np.seterr(divide='ignore', invalid='ignore') #to ignorewarning using missing value for calculation e.g (1.0 or 1/NA)
var_daily = (svar * mask).sum(axis=0) / mask.sum(axis=0)
np.seterr(divide='warn', invalid='warn') #to reactivate warning using missing value for calculation
idx = np.where(np.isnan(var_daily) == True)
var_daily[idx] = spval
var_daily = np.ma.masked_values(var_daily, spval)
nc.variables[outvarname][nday,:,:] = var_daily
nday = nday + 1
# dataset.close()
nc.close()
Nice to meet you, let's keep in touch!