## l3mapgen georeferencing

Use this Forum to find information on, or ask a question about, NASA Earth Science data.
OB.DAAC - SeanBailey
User Services
Posts: 1429
Joined: Wed Sep 18, 2019 6:15 pm America/New_York
Been thanked: 1 time

### l3mapgen georeferencing

> 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

Tags:

loicdtx
Posts: 9
Joined: Wed Nov 30, 2016 8:19 am America/New_York

### l3mapgen georeferencing

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
Subject Matter Expert
Posts: 142
Joined: Tue Feb 09, 2021 8:19 am America/New_York

### l3mapgen georeferencing

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/33r10.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/33r1.76527037539  13.13061819echo "-72 33"|gmt mapproject -JA-97/18/20 -R-122/3/-72/33r20  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/33r21.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.psecho -e "-122 3\n-72 3"|gmt psxy -JA -R -P -Ap -W2,255/0/0 -O -K>>map.psecho -e "-122 33\n-72 33"|gmt psxy -JA -R -P -Ap -W2,255/0/0 -O -K>>map.psecho -e "-122 3\n-122 33"|gmt psxy -JA -R -P -Am -W2,255/0/0 -O -K>>map.psecho -e "-72 3\n-72 33"|gmt psxy -JA -R -P -Am -W2,255/0/0 -O>>map.psgs -sDEVICE=pngalpha -sOutputFile=- -q -r100 -g2200x1461 -<map.ps >map.png`
I hope someone finds this helpful.

Regards,
Norman

loicdtx
Posts: 9
Joined: Wed Nov 30, 2016 8:19 am America/New_York
`import pyprojfrom affine import Affineimport rasteriofrom rasterio.crs import CRSimport netCDF4 as ncimport numpy as npinput_filename = 'T2014027.L3m_DAY_CHL.nc'proj4string = '+proj=laea +lat_0=18 +lon_0=-97'# Open netcdf connectioncon = nc.Dataset(input_filename)# Array shape and resolutionres = con.geospatial_lon_resolution * 1000height = con.number_of_lineswidth = con.number_of_columns# Get dictionary representation of proj4stringcrs = CRS.from_string(proj4string)# Define pyproj transformation objectp = 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 computedsw_ll = (con.westernmost_longitude, con.southernmost_latitude)# Southest point of the projected extent should be at lon_0 (crs definition)/southernmost_latitudesouth_max_ll = (crs['lon_0'], con.southernmost_latitude)# Convert coordinate pairs to projected CRSx_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 rasteriogeo_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 filechlor_a = con.variables['chlor_a'][:]con.close()# # Create georeferenced file and write array to itwith rasterio.open('T2014027.L3m_DAY_CHL.tif', 'w', **geo_dict) as dst:    dst.write_band(1, chlor_a.astype(rasterio.float32))`