1

I have some SMOS soil moisture data packaged in a netcdf format distributed by "cp34-bec". I have been able to use GDAL to extract this data as a raster layer, using a projection i think is appropriate and a geotransform that is present in the netcdf.

The problem is that, while the image should span an area on earths surface 5x8 degrees in size, it instead spans an area that is 5x8 meters in size. Essentially, if these units were being treated as degrees in Arcmap and Qgis instead of meters, everything would be great. I have looked through the GDAL metadata commands and haven't been able to find one that allows me to identify some kind of unit confusion.

the metadata says it is an EASE_ML grid, which can be read about here, and the projection i have assigned can be found at spatialreference.org. I have considered that the projection needs to be amended to reflect the unit preference, but haven't been able to make a meaningful change on my own. Any help would be great!

Since it was asked about, here is the whole function i use to write gdal datasets to tiff.

import gdal
from convert_dtype import *

def _gdal_dataset_to_tif(gdal_dataset, outpath, cust_projection = None,
                     cust_geotransform = None, force_custom = False,
                     nodata_value = None):
"""
This function takes a gdal dataset object as returned from the
"_extract_HDF_layer_data" OR "_extractNetCDF_layer_data functions
and writes it to tif with either the embedded projection and
geotransform or custom ones. This function should be wrapped in another
function for a specific datatype.

:param gdal_dataset:            a gdal.Dataset object
:param outpath:                 output filepath for this dataset (tif)
:param cust_projection:         a projection string, see datatype_library
:param cust_geotransform:       a geotransform array, see datatype_library
:param force_custom:            if True, forces the custom geotransform and
                                projections to be used even if valid
                                geotransforms and projections can be read
                                from the gdal.dataset. If False, custom
                                projections and geotransforms will be ignored
                                if valid variables can be pulled from the
                                gdal.dataset metadata.
:param nodata_value:            the value to set to Nodata

returns the local system filepath to output dataset.
"""

# set up the projection and geotransform
if force_custom is True:
    projection = cust_projection
    geotransform = cust_geotransform
else:

    gdal_projection = gdal_dataset.GetProjection()

    # only uses the custom projection if gdal metadata is bad
    if gdal_projection == "":
        projection = cust_projection
    else:
        projection = gdal_projection

    gdal_geotransform = gdal_dataset.GetGeoTransform()

    # only uses the custom geotransform if gdal geotransform is default (bad)
    if gdal_geotransform == (0.0, 1.0, 0.0, 0.0, 0.0, 1.0):
        geotransform = cust_geotransform
    else:
        geotransform = gdal_geotransform


# print update
print("using projection {0}".format(projection))
print("using geotransform {0}".format(geotransform))

# set up the numpy array
numpy_array = gdal_dataset.ReadAsArray()
shape = numpy_array.shape

# determine its shape
if len(shape) == 2:
    xsize = shape[1]
    ysize = shape[0]
    numbands = 1
elif len(shape) == 3:
    xsize = shape[2]
    ysize = shape[1]
    numbands = shape[0]
else:
    raise Exception("cannot write 1 dimensional data to tif")

# create the tiff
gtiff = gdal.GetDriverByName("GTiff")
outdata = gtiff.Create(outpath, xsize, ysize, numbands, convert_dtype(numpy_array.dtype))
outdata.SetProjection(projection)
outdata.SetGeoTransform(geotransform)

# write each band
for i in range(numbands):
    outraster = outdata.GetRasterBand(i+1)
    outraster.WriteArray(numpy_array, 0, 0)
    if nodata_value is not None:
        outraster.SetNoDataValue(nodata_value)
    outraster.FlushCache()

return outpath

Aside: Please forgive my small number of sources. The reputation system prevents me from including more than two hyperlinks. I have also been unable to appropriately tag this question.

8
  • You'll get more rep, which leads to more ability, it's only because it's your first question. If you're sure that the spatial reference is correct use Define Projection tool (or right click in catalog, go to properties and set the spatial reference there) you can search in the XY Coordinate System dialog with 3410 to get the right one. Jul 1, 2015 at 0:06
  • This is really a comment but likewise lack the reputation points to comment. What is the GDAL command to extract the raster layer?
    – elil
    Jul 1, 2015 at 0:35
  • What's the value of projection? My go-to is from osgeo import osr; osr.UseExceptions(); srs = osr.SpatialReference(); srs.ImportFromEPSG(3410); projection = srs.ExportToWkt()
    – Mike T
    Jul 1, 2015 at 0:45
  • 1
    Also the units for EPSG:3410 is metres, not degrees.
    – Mike T
    Jul 1, 2015 at 0:46
  • The value of projection is the string read from a .prj file, but I just checked and discovered the .prj file that is distributed by spatialreference.org is NOT identical to the string returned by your code method. However, ive just tried it out using the osr library and the exact same output results.
    – Jwely
    Jul 1, 2015 at 0:59

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.