Page History
| Lucidchart | ||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
GDAL can access OPeNDAP. However, it can handle dataset that has 2 dimensions only. Use netCDF or pydap for generic access.
dimension. If you want to specify constraint expression, it should have "x" and "y" as explained in the GDAL user guide. The following code is from frmts/dods/dodsdataset2.cpp:
| Code Block | ||
|---|---|---|
| ||
/* -------------------------------------------------------------------- */
/* For now we hard code to assume that the two dimensions are */
/* ysize and xsize. */
/* -------------------------------------------------------------------- */
if( poArray->dimensions() < 2 )
{
throw Error("Variable does not have even 2 dimensions. For now this is req\
uired." );
}
int nXDir = 1;
int nYDir = 1;
int iXDim = GetDimension( oCE, "x", &nXDir );
int iYDim = GetDimension( oCE, "y", &nYDir );
if( iXDim == -1 || iYDim == -1 )
{
throw Error("Missing [x] or [y] in constraint." );
} |
Therefore, you need to specify subsetting carefully. Here's a working Python example.
| Code Block | ||
|---|---|---|
| ||
from osgeo import gdal
dods_dr = gdal.GetDriverByName('DODS')
if dods_dr is None:
print('DODS driver is missing.')
ds = gdal.Open('https://opendap.larc.nasa.gov/opendap/hyrax/CERES/SYN1deg-1Hour/Terra-Aqua-MODIS_Edition4A/2018/01/CER_SYN1deg-1Hour_Terra-Aqua-MODIS_Edition4A_403405.20180101.hdf?sza[0:1:0][y][x]')
print(ds.RasterCount)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
print(arr) |
Please pay attention to the 'x' and 'y' in URL in the above code.
Here's the basic GDAL to GeoTIFF code.
| Code Block | ||
|---|---|---|
| ||
# Continue from OPeNDAP to GDAL Workflow code
arr = band.ReadAsArray()
[cols, rows] = arr.shape
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create("CER.tif", rows, cols, 1, gdal.GDT_UInt16)
geotransform = ([-180, 1, 0, 90, 0, -1 ])
outdata.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
outdata.SetProjection(srs.ExportToWkt())
outdata.GetRasterBand(1).WriteArray(arr)
# Saves to disk.
outdata.FlushCache()
outdata = None
band=None
ds=None
|
We can confirm the correctness using QGIS for GeoTIFF and Panoply for OPeNDAP.
| Panoply (longitude and latitude are sub-setted as well) | QGIS (-180 to 180) Grey scale: white is 90, black is 19. |
|---|---|
|
|
Overview
Content Tools