l3mapgen georeferencing

[OB.DAAC - SeanBailey]

> Any idea what's happening.

Yes, I think I may.
Your assumption that the northwest corner of the dataset is -122, 33 is incorrect. 
l3mapgen requires a rectangular shape when it generates it's output raster.  The
input bounding coordinates are used to define that rectangle such that the output
will completely contain the defined region of interest.  In your case, you have defined
a non-rectangular projection, so the output raster is larger to accommodate the ROI.
For this particular scene and projection, the lower left is -122, 3; lower right -72,3 as
you specified (happens that the widest part of the image matches the input).
However, the UPPER corners of the box are NOT -122, 33 and -72,33:

Sean

[loicdtx]

Yes, that makes perfect sense.
I therefore tried to georeference the array using the south-west corner coordinates but the result is off too. Which also makes sense, since for the 3 degrees parallel to be entirely contained in the bounding box, the south-west corner must be below 3 degrees.
Unfortunately I could not find any, even inelegant, way to compute the real coordinates of any of the corners.
Any idea?

Cheers,
Loïc

[OB WebDev - norman]

Well, here is one inelegant way to get the corner coordinates you seek.  I am using
GMT in my example since that is what I am most familiar with, but the same logic should
work using any other map projection software like Proj4.

# First, naively specify a Lambert azimuthal equal area projection with a
# lower left corner at 122W by 3N and an upper right corner at 72W by 33N
# with the projection centered at 97 W by 18 N.  In GMT parlance this looks
# like:
#    -R-122/3/-72/33r -JA-97/18/20
# (The final "20" just says I want my map to be 20 units wide.)

# Intuitive knowledge of this projection in the northern hemisphere lets me
# know that the lowest part of the 3N parallel will be at the central
# meridian of the projection i.e. 97W, so I project that coordinate to
# x and y in my naive map space.
echo "-97 3"|gmt mapproject -JA-97/18/20 -R-122/3/-72/33r
10.8826351877  -0.589113989394

# I also know that the highest part of the 33N parallel will be at
# whichever edge is farther from the central meridian.  In this case the
# central meridian of the projection is midway between the western and
# eastern boundaries, so I can choose either one.
echo "-122 33"|gmt mapproject -JA-97/18/20 -R-122/3/-72/33r
1.76527037539  13.13061819
echo "-72 33"|gmt mapproject -JA-97/18/20 -R-122/3/-72/33r
20  13.13061819
# As you see the y-coordinate is the same in both cases.

# Finally, the rightmost part of my desired map will be at the southeast
# corner of my region of interest.
echo "-72 3"|gmt mapproject -JA-97/18/20 -R-122/3/-72/33r
21.7652703754  0

# Now that I have the x and y boundaries of the map I want, I can do the
#  inverse (-I) projection to get the lower left and upper right
# geographical coordinates I want.
# Lower left:
echo "0 -0.589113989394"|gmt mapproject -JA-97/18/20 -R-122/3/-72/33r -I
-121.93314575  1.66297305618
# Upper right:
echo "21.7652703754 13.13061819"|gmt mapproject -JA-97/18/20 -R-122/3/-72/33r -I
-67.2236919686  32.2700636011

# The above coordinates now let me construct the actual map I desire.
# For a 20-unit-wide map, the GMT arguments are (keeping way more digits
# than needed so you can see where the numbers come from):
# -R-121.93314575/1.66297305618/-67.2236919686/32.2700636011r -JA-97/18/20

# The attached map was generated as follows with extra commands to draw
# in the region of interest in red.
gmt pscoast -JA-97/18/20 -R-121.933/1.663/-67.224/32.270r -P -X1 -Y1 -G200 -S100 -B2g1 -Di --PROJ_LENGTH_UNIT=i --FORMAT_GEO_MAP=DF -K>map.ps
echo -e "-122 3\n-72 3"|gmt psxy -JA -R -P -Ap -W2,255/0/0 -O -K>>map.ps
echo -e "-122 33\n-72 33"|gmt psxy -JA -R -P -Ap -W2,255/0/0 -O -K>>map.ps
echo -e "-122 3\n-122 33"|gmt psxy -JA -R -P -Am -W2,255/0/0 -O -K>>map.ps
echo -e "-72 3\n-72 33"|gmt psxy -JA -R -P -Am -W2,255/0/0 -O>>map.ps
gs -sDEVICE=pngalpha -sOutputFile=- -q -r100 -g2200x1461 -<map.ps >map.png

I hope someone finds this helpful.

Regards,
Norman

[loicdtx]

Not inelegant at all :) For some reasons I could not visualize that -122,3 would give me the western edge of the projected extent... So thanks for showing your reasoning, it helped a lot even though gmt isn't part of my workflow.

For future reference, I'm posting below the python way to get a referenced geotiff from the l3mapgen netcdf file.

Cheers,
Loïc

import pyproj
from affine import Affine
import rasterio
from rasterio.crs import CRS
import netCDF4 as nc
import numpy as np

input_filename = 'T2014027.L3m_DAY_CHL.nc'
proj4string = '+proj=laea +lat_0=18 +lon_0=-97'

# Open netcdf connection
con = nc.Dataset(input_filename)

# Array shape and resolution
res = con.geospatial_lon_resolution * 1000
height = con.number_of_lines
width = con.number_of_columns

# Get dictionary representation of proj4string
crs = CRS.from_string(proj4string)

# Define pyproj transformation object
p = pyproj.Proj(proj4string)

# Two pairs of coordinates are required to compute corner coordinates in
# projected CRS
# South-west corner of the unprojected extent from which 'xmin'
# will be computed
sw_ll = (con.westernmost_longitude, con.southernmost_latitude)
# Southest point of the projected extent should be at lon_0 (crs definition)/southernmost_latitude
south_max_ll = (crs['lon_0'], con.southernmost_latitude)

# Convert coordinate pairs to projected CRS
x_min, y_dummy = p(*sw_ll)
x_dummy, y_min = p(*south_max_ll)
y_max = y_min + height * res

# Build metadata dictionary required by rasterio
geo_dict = {'affine': Affine(res, 0.0, x_min,
                             0.0, -res, y_max),
            'height': height,
            'width': width,
            'crs': crs}

geo_dict.update(driver = u'GTiff', dtype = rasterio.float32, count = 1, nodata = -32767)

# # Read an array of the netcdf file
chlor_a = con.variables['chlor_a'][:]
con.close()

# # Create georeferenced file and write array to it
with rasterio.open('T2014027.L3m_DAY_CHL.tif', 'w', **geo_dict) as dst:
    dst.write_band(1, chlor_a.astype(rasterio.float32))