1

I have an array made from MODIS Terra Daily snow data that I would like to export for use in ArcMap. How would I output the 'probablility_of_snow array' array near the bottom of this post to an image to use in ArcMap?

Here are some of the steps I followed to compute the array:

 from pyhdf.SD import SD, SDC
 import numpy as np
 import matplotlib.pyplot as plt
 import gdal

 #Get 2001 data
 file = SD(MOD10A1_A2001001_h09v04_006_2016091183527_HEGOUT.hdf], SDC.READ)
 datasets_dic = file.datasets()
 print ('datasets_dictonary\n', datasets_dic)
 print ('file info', file.info())

datasets_dictonary:

 {'NDSI_Snow_Cover': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 21, 0), 'NDSI_Snow_Cover_Basic_QA': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 21, 1), 'NDSI_Snow_Cover_Algorithm_Flags_QA': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 21, 2), 'NDSI': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 22, 3), 'Snow_Albedo_Daily_Tile': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 21, 4), 'orbit_pnt': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 20, 5), 'granule_pnt': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 21, 6)}

File.info

 file info (7, 5)

I get the data for several years and I assign each year's data to an array:

 sds_obj = file.select('NDSI_Snow_Cover') 
 Array_2001 = np.array(sds_obj.get()) 
 print (Array_2001)

The data looks like this for each year:

 [[ 84 85 80 ..., 250 250 250]
 [ 82 82 85 ..., 250 250 250]
 [ 78 82 82 ..., 250 250 250]
 ..., 
 [250 250 10 ..., 250 250 250]
 [250 250 250 ..., 250 250 250]
 [250 250 250 ..., 250 250 250]

I do some data cleaning and math (not all shown here):

 probablility _of_snow = np.divide(snow_days,sample_size)

This is the array I want to export to make a map in ArcMap:

 probablility_of_snow array  :

  [[ 0.75        0.75        0.71428571 ...,  0.71428571  0.71428571
  0.71428571]
  [ 0.625       0.625       0.625      ...,  0.71428571  0.71428571
  0.71428571]
  ..., 
  [ 0.42857143  0.5         0.5        ...,  0.71428571  0.71428571
  0.71428571]
  [ 0.42857143  0.5         0.375      ...,  0.71428571  0.71428571
  0.71428571]]

Here is some info from the MODIS file:

 #Get transformation and projection info
 width = 1038
 height = 527
 file = SD(MOD10A1_A2001001_h09v04_006_2016091183527_HEGOUT.hdf, SDC.READ)
 sds = gdal.Open(file, gdal.GA_ReadOnly).GetSubDatasets()
 geoInf = gdal.Open(sds[0][0])
 geoT = geoInf.GetGeoTransform()
 proj = geoInf.GetProjection()


 print(geoT)

 (-118.1513669444455, 0.0045054006101485754, 0.0, 44.71875, 0.0, -0.004505401117436153)


 print (proj)

 GEOGCS["Unknown datum based upon the Clarke 1866 ellipsoid",DATUM["Not specified (based on Clarke 1866 spheroid)",SPHEROID["Clarke 1866",6378206.4,294.9786982138982,AUTHORITY["EPSG","7008"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]

I don't understand what to do next to create and export an image to use in ArcMap.

  • If you also want to ask about Google Earth Engine then feel free to do that in a separate question. – PolyGeo Mar 9 '18 at 21:40
1

Assuming you have an ArcGIS license you can use the NumpyArrayToRaster function: https://pro.arcgis.com/es/pro-app/arcpy/functions/numpyarraytoraster-function.htm


In case you want to use gdal it is going to take a few more lines. Here is an example:

Note: gdal needs you to specify a Geotransform and a Projection, as well as the number of columns and rows of the dataset, when creating a new raster. I assume here they are your geoT and proj variables, respectively, and that the size of the raster is the same one as the geoInf data set.

# I assume you already imported all the necessary modules
# Your processing goes here

out_fn = r"C:\...\output.tif"
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(out_fn, geoInf.RasterXSize, geoInf.RasterYSize, 1,
                       gdal.GDT_Float32)
out_ds.SetProjection(proj)
out_ds.SetGeoTransform(geoT)

out_band = out_ds.GetRasterBand(1)
out_band.SetNoDataValue(-99)

out_band.WriteArray(probablility _of_snow)
out_band.FlushCache()
out_band.ComputeStatistics(False)

del geoInf, out_ds

Further note

You can read datasets (including .hdf files) as NumPy arrays using gdal. For example:

ds = gdal.Open(r"C:\...\input.tif")
array = ds.GetRasterBand(n).ReadAsArray()  # 2d array

ds = gdal.Open(r"C:\...\input.tif")
array = ds.ReadAsArray()  # 3d array
  • Is there another way to do it using something like gdal? I would like to work on my project at home where I do not have access to ArcGis. – Steve Obert Mar 10 '18 at 2:03
  • @SteveObert I'll update my question then – Marcelo Villa Mar 10 '18 at 3:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.