How can I export a Numpy array made from MODIS data for use in ArcMap?

Question

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.

• PolyGeo: If you also want to ask about Google Earth Engine then feel free to do that in a separate question.

Answer 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

• Steve Obert: 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.
• Marcelo Villa: @SteveObert I'll update my question then