Numpy masked array to raster

Question

I am trying to transform a masked array, result of a NETCFD4 file, to a raster.

The array is saved as a masked array, with shape (500,805), and the mask is the value 32768.0. The array looks like:

import os
import gdal, ogr, osr
import numpy as np
import fiona as fn
import pandas as pd
import geopandas as gpd
import rasterio as rt
import scipy as sp
import matplotlib.pyplot as plt


masked_array(
 data=[[--, --, --, ..., 16.31, 16.39, 16.35],
       [--, --, --, ..., 16.31, 16.42, 16.34],
       [--, --, --, ..., 16.29, 16.330000000000002, 16.39],
       ...,
       [--, --, --, ..., 27.650000000000002, 28.54, 28.57],
       [--, --, --, ..., 26.64, 28.59, 27.92],
       [--, --, --, ..., 27.84, 27.830000000000002, 27.54]],
 mask=[[ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       ...,
       [ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False]],
 fill_value=32768.0)

I have written the following method, based on many posts from around. However, it fails whenever I save the file as a raster. The maximum value I get, is the one for the fill = 32768.0

def array2raster(newRasterfn, rasterOrigin, pixelWidth, pixelHeight, array, path):

   reversed_arr = array[::-1]    
   cols = reversed_arr.shape[1]
   rows = reversed_arr.shape[0]
   originX = rasterOrigin[0]
   originY = rasterOrigin[1]

   driver = gdal.GetDriverByName('GTiff')
   outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float64)
   outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
   outband = outRaster.GetRasterBand(1)                                                                                                                              

   outband.WriteArray(reversed_arr)
   outband.SetNoDataValue(-999)
   outRasterSRS = osr.SpatialReference()
   outRasterSRS.ImportFromEPSG(4326) # http://spatialreference.org/ref/epsg/32614/
   outRaster.SetProjection(outRasterSRS.ExportToWkt())
   outband.FlushCache()

I call the function like this.

jan_1950raster = array2raster('jan1950' + '.tif', [lon,lat], cell_size, cell_size, jan_1950, gis)

How can I proceed from here?

I have changed outband.SetNoDataValue(-32768.0) but the problem remains, I have a raster with these values.

Answer 1

In your method, change it like the following lines and that will solve the problem:

def array2raster(newRasterfn, rasterOrigin, pixelWidth, pixelHeight, array, path):
    """ This function works fine, until the raster mask. It does the work whatsoever. 
        Not to be used unless necessary. """

    # Create a local copy of the array
    temp = array.copy()

    # Given that the array is a masked array, the mask is "non existent".
    # Consequently, all values bigger than the maximum value (temp > np.max(temp))
    # Can be stored as "nan".
    temp[np.where(temp > np.max(temp))] = np.nan 

In the process of writing the raster, outband.SetNoDataValue will recognize the nan as the data to set to "whatever the number" you want.

Answer 2

As @Daniel mentioned, you have to specify a value for the masked positions in your masked array. A straightforward way of doing this is calling the .filled() method on said masked array. In your case you could write:

# assume arr is the masked array
arr = arr.filled(-999)  # you can pass an arbitrary scalar to this method.

Answer 3

I am working on a similar problem (wanting to work with rasterized netcdf data) and I think you may have made the same wrong assumption as me, but without information on what instrument your data comes from I can't confirm that.

In my case, I assumed the data was already in some kind of grid projection because it came in rows and columns. Instead the 'pixel grid' (really more of a polygon mosaic) is rotated 45 degrees and different passes of the same area don't overlap. This may not be the case for your data so you should doublecheck.

I just wanted to point out that in your example (even though I don't really follow along 100%) you seem to set the origin lat and lon in some way that isn't explained in the example. You should probably instead use the georeference data available in the nc file. In my file it's located in :

/PRODUCT/SUPPORT_DATA/GEOLOCATIONS/latitude_bounds

which is accessed using :

ncfileobject['PRODUCT']['SUPPORT_DATA']['GEOLOCATIONS']['latitude_bounds'][:]

(and the longitude_bounds are on the same level)

For my own issue, I am still working on a solution. Please let us know what the instrument name is, so we can help further.