# X and Y after Albers lat and lon transformation

that is probably a very obvious question but it is melting my brain:

when transforming lat and lon to [x,y] (for instance using this code https://gist.github.com/RandomEtc/476238) I get the following: { x: -0.27, y: -1.12 }

What are x and y? Are these radians? How can I use x and y to find the exact point for that lat and lon in may 2D projected map?

Edit:

This is the map I am working with: http://www.jetplan.com/weather/data/maps/eutball06.gif">

To use any formula, you need the parameters of the Albers Equal Area projection:

latitude and longitude of origin first and second parallel

You can guess the longitude from the only meridian that is straightly vertical, but the others are more difficult to get.

You might be better off by using a georeferencing tool. The georeferencer inside QGIS does a good job.

### Project Latitude-Longitude Coordinates to x - y

Project latitude-longitude coordinates to x - y coordinates by specifying a map projection. Then, display the projected coordinates on a map.

To do this, first specify the latitude and longitude coordinates of landmarks in Boston. Specify the coordinates in the NAD83 geographic CRS.

Then, import a GeoTIFF image of Boston as an array and a map cells reference object. Get information about the map projection by querying the ProjectedCRS property of the reference object. Verify that the geographic CRS underlying the projected CRS is NAD83.

Project the latitude-longitude coordinates to x - y coordinates using the same projected CRS as the GeoTIFF image.

Display the GeoTIFF image and the projected coordinates on the same map. Change the marker symbol and color of the coordinates so they are more visible. Then, add axis labels.

## GeographicLib

Comment has been marked as spam.
View and moderate all "Help" comments posted by this user

I have a system which needs to display events on a map, the map extents are small enough that I can assume a Cartesian model (affine transformation) for the points. But I need the map georeferenced, and my users want to be able to load a GeoTIFF and not manually provide map corner coordinates.

GDAL is able to do this directly from a TIFF file, by reading the gdalinfo source code I've been able to directly access the map corner coordinates that are printed by gdalinfo. But GDAL is also a huge footprint with many dependencies and merely loading its library opens TCP server ports, which is not acceptable for my project. So I'd much prefer to use GeographicLib for conversion of the georeference information.

I'm able to automatically read the GeoTIFF tag structure, create instances of the projection classes, and use the Reverse function to get latitude and longitude. For TransverseMercator this is working very well, but for other schemes the results are significantly different from what GDAL indicates. I'd appreciate an expert set of eyes to help me find the discrepancy.

The calls I'm making to GeographicLib for this file are as follows (C#, substituted actual values of the parameter variables):

the constructed projection is

As you see, these vary from the GDAL values by about 7 arc-minutes latitude and 12 arc-minutes longitude.

Could I be using the wrong k (radial scale)? Is it because the datum [EPSG 6267] isn't accounted for?

## Compatibility Considerations

### Minvtran will be removed

Not recommended starting in R2020b

The minvtran function will be removed in a future release. In most cases, use the projinv function instead. If the mapprojection field of the current map axes or specified map projection structure is 'globe' , then use the ecef2geodetic function instead.

This table shows some typical uses of the minvtran function and how to update your code to use the projinv function instead.

This table shows some typical uses of the minvtran function and how to update your code to use the ecef2geodetic function instead.

If the value of the mapprojection field of the map projection structure is listed in this table, then you must specify all linear units in meters, including coordinates and map projection structure fields. After you have projected the coordinates, you can convert them to other units using the unitsratio function.

[x,y] = mfwdtran(lat,lon)
[x,y,z] = mfwdtran(lat,lon,alt)
[. ] = mfwdtran(mstruct. )

[x,y] = mfwdtran(lat,lon) applies the forward transformation defined by the map projection in the current map axes. You can use this function to convert point locations and line and polygon vertices given in latitudes and longitudes to a planar, projected map coordinate system.

[x,y,z] = mfwdtran(lat,lon,alt) applies the forward projection to 3-D input, resulting in 3-D output. If the input alt is empty or omitted, then alt = 0 is assumed.

[. ] = mfwdtran(mstruct. ) requires a valid map projection structure as the first argument. In this case, no map axes is needed.

[x,y] = mfwdtran(lat,lon)
[x,y,z] = mfwdtran(lat,lon,alt)
[. ] = mfwdtran(mstruct. )

[x,y] = mfwdtran(lat,lon) applies the forward transformation defined by the map projection in the current map axes. You can use this function to convert point locations and line and polygon vertices given in latitudes and longitudes to a planar, projected map coordinate system.

[x,y,z] = mfwdtran(lat,lon,alt) applies the forward projection to 3-D input, resulting in 3-D output. If the input alt is empty or omitted, then alt = 0 is assumed.

[. ] = mfwdtran(mstruct. ) requires a valid map projection structure as the first argument. In this case, no map axes is needed.

## Overlaying Geographic Data on GeoTIFF Images

QUESTION: I have am trying to display an AVHRR NDVI data image in GeoTIFF format on a map projection in IDL. No matter what I do to georegister the image, the continental outlines are just slightly off. According to the published documentation, this image is in an Albers Conical Equal Area projection, with the corner pixels (from lower-left and in clockwise fashion) located at these latitude/longitude locations:

You can see from the image below that the continental outlines I put around the African continent are almost correct, but not exactly correct. The error is especially noticable around Madegascar, but you can see it along the west coast of Africa, as well. Can you help me straighten this out?

 This image is almost registered with the map projection, but not quite.

ANSWER: Because this is a GeoTIFF image, the image itself contains information that is suppose to allow the end-user to associate a map projection data coordinate space with the image itself. The data, which is called metadata, is obtained from the GeoTIFF image at the time the image is read into IDL by supplying a GEOTIFF keyword to the READ_TIFF function that reads the image data. A structure, containing GeoTIFF tags is returned by the keyword.

So you can play along if you like, I have made a GeoTIFF image available to you on my web page. Extract the GeoTIFF image to the same directory in which you will run the example IDL code, which I have named tiffoverlay.pro.

To open the GeoTIFF file and read it into IDL, type these commands.

Notice we have collected the GeoTIFF metadata into an IDL structure variable named geotag, and that we have reversed the Y direction of the image. The latter is almost always required when working with TIFF images in IDL, since the TIFF standard sets the (0,0) point in the upper-left corner of the image, and IDL almost always prefers to use the lower-left corner of the image for that purpose. I always reverse the actual image data, rather than just display the image "upside down" by setting the !Order system variable, since I often want to interact with images, and I want to avoid the confusion that is inevitable when my “location in the image” is different from my “display of the image” in the graphics window.

The GeoTIFF meta data is available as geotags, which are represented in IDL as fields in the geotag structure.

Unfortunately, most of this information doesn't mean much to you, probably, and the IDL documentation doesn't help you sort it out. To understand it, you are probably going to have to go to somewhere like RemoteSensing.org, where you can find documentation on the GeoTIFF format itself. The piece of software I have found most helpful in understanding GeoTIFF files is listgeo, which provides a dump of the GeoTIFF metadata in a human readable form. On my Windows machines, I use a GUI wrapper to listgeo, named listgeo_GUI. It is not elegant software, but I do find it essential.

To understand the specific geotags you find in your GeoTIFF file, you are going to have to search through the GeoTIFF file specifications. I usually search for the specific geotag name in this document. The information I need is typically found in the Appendices of this document.

Here is a slightly edited portion of the listgeo printout for this file. It will help you see how you can find the same information in the GeoTIFF tags directly from within IDL.

We see that the image is in an Albers Equal Area map projection. (We can obtain this information from the metadata by translating the PROJECTIONGEOKEY and PROJCOORDTRANSGEOKEY geotags--or fields in the geotag structure--according to the information in the GeoTIFF specifications.) Other information essential for setting up the map projection is also avaiable. For example, the Albers Equal Area projection requires two standard parallels and these are provided in the PROJSTDPARALLEL1GEOKEY and PROJSTDPARALLEL2GEOKEY geotags. The datum used for this projection is WGS84 (a translation of the GEOGRAPHICTYPEGEOKEY geotag, 4326). We can find the center of the map projection in the PROJNATORIGINLATGEOKEY and PROJNATORIGINLONGGEOKEY geotags.

In addition to knowing what particular map projection was used for the image, we have to know where, exactly, the image is located in that map projection. This information is also available in the geotags, but you have to know how to parse it. Of particular importance are the tie points (MODELTIEPOINTTAG) and the image scale (MODELPIXELSCALETAG). The tie points usually (but not always) locate the upper left corner of the image in Cartesian coordinates. This is simply a flat 2D surface upon which we cast the map projection and upon which we draw a grid in arbitrary units. In this case, the units are meters (from the PROJLINEARUNITSGEOKEY tag) and the extent of the grid is the diameter of the Earth. This arbitrary coordinate system is sometimes called the UV coordinate system, where U is the X dimension (often associated with longitude) and V is the Y dimension (often associated with latitude).

It is in this coordinate system that the corners and center of the image are reported in the listgeo information above. But we can calculate the corner values ourselves from within IDL. Here is the code to do it. In this code, you will see that I have moved the origin of the image from the upper-left to the lower-left. This is not strictly necessary, but I do it for my own sanity. I have been working with IDL so long that I just assume the origin is at the lower-left and I make too many mistakes if it is not. The end points are found by using the tie-points as the starting coordinates and then just multipling the size of the image in each dimension by the scaling factor for that dimension.

You can see from the print-out, that this calculation is identical to the corners reported by listgeo.

The next step is to find out where the corners of the image are in latitude and longitude space, rather than UV space. To put our image into the UV space we had to transform the data via a map projection. (In other words, our image represented some portion of a 3D surface, and we transformed that, via a map projection, into the 2D surface of the UV coordinate system.) To get back into a 3D coordinate system, represented by latitude and longitude coordinates, we have to perform an inverse transformation. We do so with MAP_PROJ_INVERSE, but to use it we have to establish the map projection space with MAP_PROJ_INIT. The code looks like this. (I have ordered the arguments so the coordinates start in the lower-left corner of the image and proceed clockwise, so the values can be compared with the published corner values.)

Printing these latitude/longitude values out, you can see they are not the same as the published values. Having spent a fair amount of time working with satellite image documentation, I would not be surprised to find the published results just simply wrong. It happens all the time. But, I think in this case, the published results describe the locations of the center of the corner pixels, and in IDL we need to work with the pixel edges. So, if you wanted to use the published corner values, you would have to move their locations by half a pixel in both the U and V dimensions to get an accurate map overlay.

To overlay data on the map projection, the extent of the map projection has to be exactly coincident with the image. Thus, we have to know the limits of the map projection. We can try to calculate those from the corner pixels. Choosing two opposite corners that have the maximum separation, we write code like this.

Unfortunately, you can see (from the calcuated longitude and latitude values above) that the limit doesn't exactly describe the location of the image on the map. We really need an eight-element limit array to describe the four corners of the image. Something like that exists in IDL for the MAP_SET command. (Although you must specify the left, top, right, and bottom location of the image, and not the corners.) But if we pass an eight-element array to the LIMIT keyword in MAP_PROJ_INIT it doesn't complain, but it completely ignores the values!

Let's see what happens if we continue to use the limit as we calculate above.We use MAP_PROJ_INIT, but with the LIMIT keyword set.

The MAP_PROJ_INIT function returns a map structure variable. One of the fields in this structure variable is UV_BOX. This field describes the extent of the map in UV coordinates. You will have to use this UV_BOX to set up a data coordinate system that you can draw on top of with MAP_CONTINENTS and MAP_GRID.

First, we set up the display and display the image. You will need programs from the Coyote Library if you wish to execute this part of the code.

Next we set up the UV data coordinate system. This is where we use the UV_BOX of the map structure. Notice the PLOT command is constructed in such a way that nothing is actually drawn on the display. Its purpose is simply to set the appropriate system variables that will establish the UV data coordinate system.

Finally, we can draw political boundaries and grid lines. Note that we have to pass the map structure to MAP_CONTINENTS and MAP_GRID so they have the means to convert their latitude/longitude data into the proper UV coordinates.

You see the result in the figure below. Notice that the boundaries are just ever so slightly off. You will notice this particularly around Madagascar, in the lower-right corner of the figure.

 The overlays don't quite match up. This is because we weren't able to specify the location of the image in the map projection exactly. This is a limitation of the LIMIT keyword of MAP_PROJ_INIT. It needs to be able to accept an eight-element LIMIT keyword instead of just a four-element LIMIT keyword.

Just for completeness, let's calculate an eight-element LIMIT keyword, as we could do if we were specifying the lat/lon data coordinate space with a MAP_SET command. The code would look like this, in which I am calculating the left, top, right, and bottom sides of the image.

You can see in the results below that the MAP_PROJ_INIT command takes the eight-element LIMIT, but ignores its values.

 The MAP_PROJ_INIT command takes an eight-element LIMIT, but completely ignores its values. The UV_BOX produced is completely wrong.

To reiterate, the problem is that we cannot extract the exact limit of the map projection for this image with a four-element LIMIT array. But, in fact, we already have the proper limits calculated from before. These are the coordinates that we put into variables xOrigin, xEnd, yOrigin, and yEnd. We simply use these variables to set up our data coordinate space. The code looks like this.

You see the correct results in the figure below.

 This image uses the calculated limits from the GeoTiff file to set up the data coordinate space and the result is perfect.

An alternative way to get perfect results, is to use the UV_BOX from MAP_PROJ_IMAGE via its UVRANGE keyword, rather than from the MAP_PROJ_INIT map structure. (After a great deal of effort, I learned this is how iMap manages to get the outlines around its maps perfect.) The MAP_PROJ_IMAGE command, unlike the MAP_PROJ_INIT command, does use an eight-element positioning keyword to calculate the UV_BOX. The commands are similar to those above. And the result looks identical to the figure directly above.

The UV_BOX values for the four methods are shown below. Notice that the UV_BOX coordinates are identical for the last two methods.

MAP_PROJ_INIT UV_BOX: -4717397.8 -4644644.7 4733971.2 4600353.9 8-Element UV_BOX: -19170318. -6842903.5 18106165. 7207980.3 Calculated UV_BOX: -4599918.3 -4615852.5 4616081.7 4600147.5 MAP_PROJ_IMAGE UV_BOX: -4599918.3 -4615852.5 4616081.7 4600147.5

Many thanks to Matt Savoie for his help in solving this long standing mystery.

## Coyote Graphics Map Object

The Coyote Graphics program cgGeoMap can do all of this georegistration for you and create a map object (cgMap) that can set up the map data coordinate system for you automatically. It can even read and return the image to you from the TIFF file. The commands to do so will look like this.

filename = 'AF03sep15b.n16-VIg.tif' cgDisplay, 500, 500, WID=5, Title='Outline cgMap Object', $XPOS=50, YPOS=50 alberMap = cgGeoMap(filename, IMAGE=tiffImage, /OnImage) scaled = BytScl(tiffImage, Top=253)+1 index = Where(scaled EQ 1) scaled[index] = 0B TVLCT, cgColor('ivory', /Triple), 0 cgLoadCT, 33, NColors=253, Bottom=1 cgImage, scaled, POSITION=pos, /KEEP_ASPECT Set up a map coordinate system. alberMap -> Draw Draw political boundaries and coasts. cgMap_Continents, /HIRES, Color='Dark Slate Blue',$ Map_Structure=alberMap, /COUNTRIES, /COASTS Add gridlines cgMap_Grid, /LABEL, Color='Slate Blue', Map_Structure=alberMap

Written 14 March 2008
Last Updated 22 November 2011

## Examples

Before using minvtran , it is necessary to create a map projection structure. You can do this with axesm or the defaultm function:

The following latitude and longitude data for the District of Columbia is obtained from the usastatelo shapefile:

This data can be projected into Cartesian coordinates of the Mercator projection using the projfwd function:

To transform the projected x-y data back into the unprojected geographic system, use the minvtran function:

## 9.1 Data preparation in R

2. Set working directory to the extracted folder in R under File - Change dir.
3. First we need to load the packages needed for the exercise

library(sp)
library(lattice)
library(rgeos)#gIntersection
library(raster)#to use "raster" function
library(zoo)

str(snowy)

#Clean up by deleting extraneous columns if needed
snowy <- snowy[c(-20:-38, -40:-64)]
snowy$Status <- snowy$E_schneide

#Make a spatial data frame of locations after removing outliers
coords<-data.frame(x = snowy$X_Coordina, y = snowy$Y_Coordina)
utm.crs<-"+proj=utm +zone=13 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
utm.spdf <- SpatialPointsDataFrame(coords= coords, data = snowy, proj4string = CRS(utm.crs))

dem <-raster("snowydem")
image(dem)
class(dem)
proj4string(dem)

#Now transform projections all to match DEM (i.e., Albers)
Albers.crs <-CRS("+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
snowy.spdf <-spTransform(utm.spdf, CRS=Albers.crs)

sublette.df <- as.data.frame(sublette.spdf)
str(sublette.df)
minx <- (min(sublette.df$x)-10830) maxx <- (max(sublette.df$x)+10830)
miny <- (min(sublette.df$y)-10830) maxy <- (max(sublette.df$y)+10830)

## create vectors of the x and y points
x <- seq(from = minx, to = maxx, by = 3610)
y <- seq(from = miny, to = maxy, by = 3610)

#Alternate bbox code for spatial points
# min max
#x -854784.4 -724665.0
#y 156859.0 247343.2

#Create vectors of the x and y points
#x <- seq(from = -854784.4, to = -724665.0, by = 3610)
#y <- seq(from = 156859.0, to = 247343.2, by = 3610)

#Create a grid of all pairs of coordinates (as a data.frame)
xy <- expand.grid(x = x, y = y)
class(xy)
str(xy)

#Identifiy projection before creating Spatial Points Data Frame
Albers.crs2 <-"+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

#NOTE: Albers.crs2 is needed because SPDF needs different projection command than spTransform above

grid.pts<-SpatialPointsDataFrame(coords= xy, data=xy, proj4string = CRS(Albers.crs2))
proj4string(grid.pts)
plot(grid.pts)
gridded(grid.pts)
class(grid.pts)

#Need to define points as a grid to convert to Spatial Polygons below
gridded(grid.pts) <- TRUE
gridded(grid.pts)
str(grid.pts)
plot(grid.pts)

#Convert grid points to Spatial Polygons in essence converting to a shapefile
gridsp <- as(grid.pts, "SpatialPolygons")
str(gridsp)
plot(gridsp)
class(gridsp)
summary(gridsp)

grid <- SpatialPolygonsDataFrame(gridsp, data=data.frame(id=row.names(gridsp), row.names=row.names(gridsp)))
class(grid)
plot(grid)
names.grd<-sapply([email protected], function(x) slot(x,"ID"))
text(coordinates(grid), labels=sapply(slot(grid, "polygons"), function(i) slot(i, "ID")), cex=0.3)

#Let's check to see if all grid cells are the same size?
summary(grid)
getSlots(class(grid))
class(slot(grid, "polygons")[[1]])
getSlots(class(slot(grid, "polygons")[[1]]))

#Check area of each cell in the sampling grid in square meters
sapply(slot(grid, "polygons"), function(x) slot(x,"area"))
#[1] 13032100

#Grid cell size converted strto square kilometers
13032100/1000000
#[1] 13.0321 is grid cell size in square kilometers

#The layer below is a mule deer HSI raster layer without disturbance from Sawyer et al. 2009 #incorporated into a layer because mule deer are considered #host for the parasite we #are investigating

nodis <-raster("snowynodis")
nodis
plot(nodis)
summary(nodis)
#Need to remove NoData from mule deer HSI layer
nodis[is.na(nodis[])] <- 0

slope = terrain(dem,opt='slope', unit='degrees')
aspect = terrain(dem,opt='aspect', unit='degrees')
dem #Now let's see metadata for each layer
slope
aspect

nlcdall <- raster("nlcd_snowy")
nlcdall #Look at raster values for the habitat layer
#Values range from 11 to 95

#Or plot to visualize categories in legend
plot(nlcdall)

#Reclassify the values into 7 groups (i.e., ll values between 0 and 20 equal 1, etc.)
m <- c(0, 19, 1, 20, 39, 2, 40, 50, 3, 51,68, 4, 69,79, 5, 80, 88, 6, 89, 99, 7)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
rc <- reclassify(nlcdall, rclmat)
plot(rc) #Now only 7 categories
rc #Now only 7 categories
class(rc)

#Check to be sure all raster have same extent for Stack creation
compareRaster(dem,slope,aspect,nodis,rc)
#[1] TRUE

#############################################################
#############################################################
##NOTE: Code in this box was simply for demonstration purposes to reduce overall time for processing during class.
##Skip this section of code if using your own data and your computer has the appropriate
## processing capabilities.

#First we will clip out the raster layers by zooming into only a few locations
plot(rc)
points(snowy.spdf)

#Code below is used to just zoom in on grid using raster layer
e <- drawExtent()
#click on top left of crop box and bottom right of crop box create zoom
newclip <- crop(rc,e)
plot(newclip)
points(snowy.spdf, col="red")

#Clip locations within extent of raster
samp_clip <- crop(snowy.spdf,newclip)
plot(newclip)
grid_clip <- crop(grid, newclip)
slope2 <- crop(slope,newclip)
aspect2 <- crop(aspect,newclip)
dem2 <- crop(dem,newclip)
HSI <- crop(nodis, newclip)

#Check to be sure all raster have same extent for Stack creation
compareRaster(dem2,slope2,aspect2,HSI,newclip)
#[1] TRUE

Your question really has nothing to do with “rectangles”, which @MPW has already pointed out don’t exist on a sphere. There is a way of looking at your problem that makes it an exercise in basic spherical trigonometry, and I’ll show you that. I do not make any claim that it’s the fastest way of getting your problem solved, though.

You have a point, let’s say $R$, on the sphere, let’s say it’s Richmond, and you want to rotate it about a fixed point $O$, let’s say that’s Orono. And your motion is a “rotation” in this sense: you have the “heading” $alpha$ from $O$ to $R$, that is the compass-direction that you start out in if you’re going in a great-circle path. And you have the distance $d$ from $O$ to $R$. Your rotation asks for the location that you’ll get to if you go a distance of $d$ with an original heading of $alpha+delta$, where $delta$ is the angle that you’re rotating things.

Now let’s draw a picture, you’ll have to do the drawing yourself. Knowing the latitude and longitude of Richmond and Orono, you see that there’s a triangle with vertex way up at the north pole $P$. The distance from $P$ to $O$ is Orono’s colatitude, that’s just the complement of Orono’s latitude. I’ll call this $c_O$. Similarly, you have the colatitude of Richmond, I’ll call this $c_R$. So you see that you have a triangle $RPO$ with legs $c_R$ and $c_O$, and up at the pole, the angle is the difference between the two cities’ longitudes. I’ll call this $lambda$, to keep Greek letters for vertex angles and Latin lower-case for lengths of sides, capital letters for points on the sphere.

So you see that you have a SAS situation, and just as in plane trigonometry your first tool to use is Law of Cosines, to get the third side of your triangle, that’s the distance from $O$ to $R$, which I’ve called $d$. Then there’s also a spherical Law of Sines for getting the vertex angle $alpha$ at $O=,$Orono. The general Laws of Cosines and Sines are $cos c=cos acos b+sin asin bcosgamma,,$ for Cosines, where $gamma$ is the angle opposite the side of length $c$. And for the Law of Sines, $frac=frac=frac,,$ where I’m sure you’ve guessed that $eta$ is to be the angle opposite the side $b$ and $alpha$ is the angle opposite the side $a$.

Now for our setup of rotating Richmond about the center Orono: We have our two legs $c_R$ and $c_O$ and our longitude-difference $lambda$, so we get $cos d = cos c_Rcos c_O + sin c_Rsin c_Ocoslambda,.$ Now, with $d$ and $lambda$ in hand, a side-and-opposite-angle pair, you can use Sines to get $alpha$. Now add your rotation-angle $delta$ to get $alpha'=alpha+delta$, and to find your Richmond