1

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.

  • 1
    Why are you specifying your nodata value as -999 instead of 32768.0? outband.SetNoDataValue doesn't change the masked values to -999, it just sets some metadata in the output raster that says treat -999 as nodata, but you haven't written any -999 values. – user2856 May 30 at 1:08
  • It would help to show the code when you call the function. – davemfish May 30 at 3:25
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.

0

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.
0

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.

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.