Как запустить opendapL4.py
Я пытаюсь загрузить спутниковые данные SST уровня 4 из PODAAC> https://opendap.jpl.nasa.gov/opendap/allData/ghrsst/data/GDS2/L4/
Я использую код python opendapL4.py, который доступен здесь https://github.com/nasa/podaac_tools_and_services/blob/master/subset_GHRSST/opendapL4.py
#! /usr/bin/env python
#
# a skeleton script to download a set of GHRSST L4 file using OPeNDAP.
#
#
# 2011.04.06 mike chin, version 0
# 2011.04.25 mike chin, version 1
# 2012.08.29 mike chin, version 2
# 2013.10.30 mike chin, version 2 (bug fix)
# 2016.01.29 mike chin, defaults to GDS2 version
##################################
# user parameters to be editted: #
##################################
# Caution: This is a Python script, and Python takes indentation seriously.
# DO NOT CHANGE INDENTATION OF ANY LINE BELOW!
# Here is the web-address (URL) for the GHRSST L4 file OPeNDAP server.
# You probably do not need to change this URL. In fact, you can
# use this URL to search for the filename of the product you want to download.
#L4URL='http://podaac-opendap.jpl.nasa.gov/opendap/allData/ghrsst/data/L4/'
#L4URL='http://opendap.jpl.nasa.gov/opendap/hyrax/allData/ghrsst/data/L4/'
#L4URL='http://seastore.jpl.nasa.gov/opendap/hyrax/allData_nc4/ghrsst/data/L4/'
#L4URL='http://opendap.jpl.nasa.gov/opendap/OceanTemperature/ghrsst/data/GDS2/L4/'
L4URL='https://opendap.jpl.nasa.gov/opendap/allData/ghrsst/data/GDS2/L4/'
# Here you set the product names in 3 parts:
# ncHead = main directory of the product.
# ncBody = main file name of the product (file name minus "yyyymmdd").
# ncTail = the part that indicates the compression method.
#
# You also specify the grid dimension of the product, which you need to
# know in advance, by setting the following 6 variables:
# nlon = x-grid (longitudes) dimension in integer.
# nlat = y-grid (latitudes) dimension in integer.
# dint = grid interval in degrees in float (assuming it's the same for x & y).
# lon0 = the smallest longitude of the grid, in degrees in float.
# lat0 = the smallest latitude of the grid, in degrees in float.
# order = how "time" ,"lon", "lat" dimensions are ordered in the L4 file,
# in array of three integers where 0=time, 1=lon, 2=lat.
#
# The code below might look complex, but it's only setting these
# 3 names (strings) and 6 variables listed above.
def strmatch(a,b):
return (a in b) and (b in a)
def parameters(product):
productlist=[ \
'mur/v4gds2', \
'mur/v4', \
'mur/v3', \
'mur/ncamerica', \
]
if strmatch( productlist[0], product.lower() ):
ncHead='GLOB/JPL/MUR/v4.1/'
ncBody='090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc'
ncTail=''
nlon=36000; nlat=17999
dint=0.01; lon0=-180.; lat0=-90.+dint
order=[0,2,1]
elif strmatch( productlist[1], product.lower() ):
ncHead='GLOB/JPL/MUR/'
ncBody='-JPL-L4UHfnd-GLOB-v01-fv04-MUR.nc'
ncTail='.bz2'
nlon=32768; nlat=16384
dint=45./(2**12); lon0=-180.+dint/2; lat0=-90.+dint/2
order=[0,2,1]
elif strmatch( productlist[2], product.lower() ):
ncHead='GLOB/JPL/MUR/'
ncBody='-JPL-L4UHfnd-GLOB-v01-fv03-MUR.nc'
ncTail='.bz2'
nlon=32768; nlat=16384
dint=45./(2**12); lon0=-180.+dint/2; lat0=-90.+dint/2
order=[0,2,1]
elif strmatch( productlist[3], product.lower() ):
ncHead='NCAMERICA/JPL/MUR/'
ncBody='-JPL-L4UHfnd-NCAMERICA-v01-fv02-MUR.nc'
ncTail='.gz'
nlon=13501; nlat=8201;
dint=0.01; lon0=-165.; lat0=-20.
order=[0,2,1]
else:
msg='\nNo such "product".\nAvailable products are:\n'
for i in range(len(productlist)):
msg=msg+' '+productlist[i]+'\n'
import sys
sys.exit(msg)
return(ncHead,ncBody,ncTail,nlon,nlat,dint,lon0,lat0,order)
# Done!!
#
# Now get out of the editor and issue the command
# opendapL4.py -h
# or
# python opendapL4.py --help
# to see how you could specify
# the dates of the files and lon-lat box from the command line.
#############################################
# You're done. No need to edit lines below #
#############################################
import sys,os
from math import floor,ceil
from optparse import OptionParser
###############
# subroutines #
###############
def ncname(body,year,yday):
(day,month)=calday(yday,year)
return('%04d%02d%02d%s'%(year,month,day,body)) # for GDS1.x and GDS2.0
def pathname(head,body,tail,y,d):
return( '%s%04d/%03d/%s%s' %(head,y,d,ncname(body,y,d),tail) )
def yearday(day,month,year):
months=[0,31,28,31,30,31,30,31,31,30,31,30,31]
if isLeap(year):
months[2]=29
for m in range(month):
day=day+months[m]
return(day)
def span(i1,i2,i3=1):
return range(i1,i2+i3/abs(i3),i3)
def isLeap(year):
flag = ( (year%4)==0) and ( not ( (year%100)==0 and (year%400)!=0 ))
return(flag)
def calday(yday,year):
months=[0,31,28,31,30,31,30,31,31,30,31,30,31]
if isLeap(year):
months[2]=29
for m in span(1,12):
months[m]=months[m]+months[m-1]
for m in span(1,12):
if (yday-1)/months[m]==0:
month=m
day=yday-months[m-1]
return (day,month)
import sys
sys.exit('ERROR calday: yearday value out of range')
def cal2mjd(year,month,day):
import math
a = (14 - month) // 12
y = year + 4800 - a
m = month + (12 * a) - 3
p = day + (((153 * m) + 2) // 5) + (365 * y)
q = (y // 4) - (y // 100) + (y // 400) - 32045
return int( math.floor(p + q - 2400000.5) )
def mjd2cal(mjd):
y=(mjd-409)//366+1860
while cal2mjd(y+1,1,1)<=mjd:
y=y+1
m=1
while cal2mjd(y,m+1,1)<=mjd:
m=m+1
d=1
while cal2mjd(y,m,d+1)<=mjd:
d=d+1
return(y,m,d)
def today(daysago=0):
import datetime
thisday=datetime.date.fromordinal(datetime.date.today().toordinal()-daysago)
return( thisday.year, thisday.month, thisday.day )
def parseoptions(L4URL):
usage = "usage: %prog [options] "
parser = OptionParser(usage)
parser.add_option("-p", "--product", help="product name [default: %default]",
dest="product", default="MUR/v4GDS2")
parser.add_option("-s", "--start",
help="start date: yyyy mm dd [default: yesterday %default]",
type="int", nargs=3, dest="date0", default=today(1))
parser.add_option("-f", "--finish",
help="finish date: yyyy mm dd [default: start date]",
type="int", nargs=3, dest="date1", default=-1)
parser.add_option("-d", "--dayspan",
help="sampling period in days [default: %default]",
type="int", dest="period", default=1)
parser.add_option("-r", "--region", "-b", "--subset",
help="limit the domain to a box given by lon_min,lon_max,lat_min,lat_max",
type="float", nargs=4, dest="box", default=(-180.,180.,-90.,90.))
parser.add_option("-g", "--gridspan",
help="sampling interval over grid [default: %default]",
type="int", dest="gridpoints", default=1)
parser.add_option("-e", "--error", help="include uncertainty error field",
action="store_true", dest="includeError",default=False)
parser.add_option("-l", "--landmask", help="include landmask field",
action="store_true", dest="includeLandMask",default=False)
parser.add_option("-i", "--ice", help="include ice concentration field",
action="store_true", dest="includeIce",default=False)
parser.add_option("-n", "--noSST", help="exclude SST field",
action="store_false", dest="includeSST",default=True)
parser.add_option("-w", "--wget", help="use wget instead of curl",
action="store_true", dest="useWget",default=False)
parser.add_option("-u", "--url", help="URL for GHRSST L4 OPeNDAP server [default: %default]", dest="URL", default=L4URL)
parser.add_option("-t", "--test", help="only prints, without issuing, the OPeNDAP commands.",
action="store_true", dest="testing",default=False)
(options, args) = parser.parse_args()
if len(args) != 0:
parser.error("incorrect number of arguments")
return( options )
# get command line options:
options=parseoptions(L4URL)
L4URL = options.URL
date0 = options.date0
if options.date1==-1:
date1 = date0
else:
date1 = options.date1
year0=date0[0]; month0=date0[1]; day0=date0[2];
year1=date1[0]; month1=date1[1]; day1=date1[2];
dday = options.period
box = list( options.box )
ig = options.gridpoints
product = options.product
p=parameters(product)
ncHead=p[0];ncBody=p[1];ncTail=p[2]
nlon=p[3];nlat=p[4];dint=p[5];lon0=p[6];lat0=p[7];order=p[8]
# find index set based on the given "box":
def boundingindex(dmin,dint,length,boundary0,boundary1):
inx0=max(int(ceil((boundary0-dmin)/dint)),0)
inx1=min(int(floor((boundary1-dmin)/dint)),length-1)
return [inx0,inx1]
[i0,i1]=boundingindex(lon0,dint,nlon,box[0],box[1])
[j0,j1]=boundingindex(lat0,dint,nlat,box[2],box[3])
if i0>i1 or j0>j1:
sys.exit('No grid point in your domain box.')
# modify the max grid indices, as necessary:
if ig>1:
i1=max(span(i0,i1,ig))
j1=max(span(j0,j1,ig))
# modify the finish day, as necessary:
if dday>1:
maxmjd=max(span(cal2mjd(year0,month0,day0),cal2mjd(year1,month1,day1),dday))
(year1,month1,day1)=mjd2cal(maxmjd)
# download size information:
print(' ')
print('First file: '+ncname(ncBody,year0,yearday(day0,month0,year0)))
print('Last file: '+ncname(ncBody,year1,yearday(day1,month1,year1)))
print(' files obtained at %d-day interval'%(dday))
print(' ')
print('Longitude range: %f to %f'%(box[0],box[1]))
print('Latitude range: %f to %f'%(box[2],box[3]))
print(' every %d pixel(s) is obtained'%(ig))
print(' ')
print('grid dimensions will be ( %d x %d )'%(len(span(i0,i1,ig)),len(span(j0,j1,ig))))
print(' ')
r=raw_input('OK to download? [yes or no]: ')
if len(r)==0 or (r[0]!='y' and r[0]!='Y'):
print('... no download')
sys.exit(0)
# form the index set for the command line:
inx=[[0,1,0],[i0,ig,i1],[j0,ig,j1]]
index=''
for i in order:
index=index+'[%d:%d:%d]'%(inx[i][0],inx[i][1],inx[i][2])
# main loop:
for mjd in span(cal2mjd(year0,month0,day0),cal2mjd(year1,month1,day1),dday):
(year,month,day)=mjd2cal(mjd)
yday=yearday(day,month,year)
cmd=L4URL+pathname(ncHead,ncBody,ncTail,year,yday)
cmd=cmd+'.nc?'
if options.includeSST:
cmd=cmd+'analysed_sst'+index+','
if options.includeError:
cmd=cmd+'analysis_error'+index+','
if options.includeLandMask:
cmd=cmd+'mask'+index+','
if options.includeIce:
cmd=cmd+'sea_ice_fraction'+index+','
cmd=cmd[0:(len(cmd)-1)] # remove the extra "," at the end.
if options.useWget:
cmd='wget "'+cmd+'" -O '+ncname(ncBody,year,yday)
else:
cmd='curl -g "'+cmd+'" -o '+ncname(ncBody,year,yday)
if options.testing: # just test command:
print(cmd)
else: # run opendap
os.system( cmd )
Я пробовал использовать терминал ubuntu и в Windows через терминал anaconda, используя команду
python opendapL4.py --start 2020 01 01 --finish 2020 12 31 --region 178 179 -19 -18
Однако я продолжаю получать сообщение об ошибке и ничего не загружается,
ошибка
Traceback (most recent call last):
File "opendapL4.py", line 282, in <module>
print('First file: '+ncname(ncBody,year0,yearday(day0,month0,year0)))
File "opendapL4.py", line 130, in ncname
(day,month)=calday(yday,year)
File "opendapL4.py", line 155, in calday
for m in span(1,12):
File "opendapL4.py", line 145, in span
return range(i1,i2+i3/abs(i3),i3)
TypeError: 'float' object cannot be interpreted as an integer
Мы очень ценим любые намеки на то, что я делаю неправильно.
С уважением,J