how to overlay psbasemap onto USGS 1:100k quads?

[Dan Scheirer]

Greetings,

I would like, but have not been able, to plot data in GMT so that I can superimpose
my PS file onto digital (PDF) versions of USGS 1:100,000 (60'x30')
topographic quadrangle maps.  The quad in question is in central
California (Cuyama), at the western edge of UTM zone 11 North,
but I think that my question would apply to many/all such
100k quad maps produced by the USGS (and perhaps the other
scale series, too).

If you have experience/interest with this, please read-on...

Thanks,
Dan Scheirer


I thought this would be easy to do because the map
legend says: "Projection and 10000 meter grid, zone 11
Universal Transverse Mercator".  I tried:
psbasemap -Ju+11/1:100000 -Ba0.5 -R-120/-119/34.5/35 -X1i -Y1i >! test.ps
It's immediately apparent that the lon-lat grid-lines are not parallel
to the edges of the paper, esp. because this is on the edge of a
UTM zone, and this makes sense to me.  If I use the -R<lower-left>/<upper-right>r
specification, I get edges to the map-box parallel to the edges of the paper, but
I do not get the full lat-lon range as specified above (thin "triangles" are missing).  
This also makes sense to me given the TM distortion, maximal near the edge of a zone.

What is unexpected to me is that the USGS quad map has the min-max longitudes
and min-max latitudes perfectly define the corners of the map, and that
these corners define edges that are parallel to the edges of the paper (parallel
to the best of my ability to measure distances on these ~archD-sized
maps).  Looking at the 10000 meter grid-lines on the quad, they are indeed
not parallel to the edges of the plotted box (nor to the edges of the paper).  
These observations suggest to me that the 100k quad maps are not actually
plotted with a UTM projection. However, when reading the other map information,
I cannot see other projection information, other than (perhaps) something to do
with aerial photograph planimetry.

The only two -J projections for geographic data that I know off the top of my head that would
make lon/lat lines be perfectly parallel to the edges of a map are: -JM|m and -JX|x...d
The equatorial Mercator option, -Jm1:100000  gives a decent overlay, but it seems funny to me
that the topo-map producers would use a projection with such distortion far away from the
Equator.  If my target plot had any more latitude range to it, I suspect that -Jm would overlay
unsuitably.

For the fake-it-as-a-linear-plot-with-the-R-coordinates-in-degrees option, -Jx1:100000d, I get
the -J syntax error message:
% psbasemap -R-120:00/-119:00/34:30/35:00 -Jx1:100000d -Ba0.5 -X1 -Y1 --PAPER_MEDIA=archD >! test.Jx.ps
psbasemap: GMT SYNTAX ERROR -J option.  Correct syntax:
....rest of the error message...

Two questions:
1) Has anyone perfectly-overlain GMT-produced *.ps files onto 1:100k
quads?  If so, with what -J projection option.  Some searching through
this list and the whole Internet didn't produce wisdom on this.  I looked
at a few 1:24,000 USGS topo quads that I have, and these seem to
have the same issue: the map legend indicates a polyconic projection
but the boundary of the maps are perfectly at the 7.5-minute lon-lat
values and these map-boundaries are perfectly parallel to the edges
of the printed map.

2) Should -Jx1:100000d produce a usage error?  -JX1d does not produce
an error, and it yields a plot with FANCY basemap type as I would
expect for a geographic-aware map.


By the way, the GeoTIFF associated with this quad has metadata (tags) indicating
that this file is indeed in UTM 11 North; when I view this in OpenEV, the lon-lat
edges of the map's outer box are indeed not parallel to the graphic window,
and the skewed-angle looks just like that of the psbasemap -Ju+11... map.

I may ask my ESRI-savvy colleagues how they would produce
an overlay map, but thought I'd try this list first.

Thanks!


     

To unsubscribe, send the message "signoff gmt-help" to listserv@...

[Mike Willis]

If you scanned the map, then the image is already projected, so you
want the linear scale. Take a look at example 28.

Cheers

On Mon, Oct 5, 2009 at 4:54 PM, Dan Scheirer <dssmailing@...> wrote:

> Greetings,
>
> I would like, but have not been able, to plot data in GMT so that I can superimpose
> my PS file onto digital (PDF) versions of USGS 1:100,000 (60'x30')
> topographic quadrangle maps.  The quad in question is in central
> California (Cuyama), at the western edge of UTM zone 11 North,
> but I think that my question would apply to many/all such
> 100k quad maps produced by the USGS (and perhaps the other
> scale series, too).
>
> If you have experience/interest with this, please read-on...
>
> Thanks,
> Dan Scheirer
>
>
> I thought this would be easy to do because the map
> legend says: "Projection and 10000 meter grid, zone 11
> Universal Transverse Mercator".  I tried:
> psbasemap -Ju+11/1:100000 -Ba0.5 -R-120/-119/34.5/35 -X1i -Y1i >! test.ps
> It's immediately apparent that the lon-lat grid-lines are not parallel
> to the edges of the paper, esp. because this is on the edge of a
> UTM zone, and this makes sense to me.  If I use the -R<lower-left>/<upper-right>r
> specification, I get edges to the map-box parallel to the edges of the paper, but
> I do not get the full lat-lon range as specified above (thin "triangles" are missing).
> This also makes sense to me given the TM distortion, maximal near the edge of a zone.
>
> What is unexpected to me is that the USGS quad map has the min-max longitudes
> and min-max latitudes perfectly define the corners of the map, and that
> these corners define edges that are parallel to the edges of the paper (parallel
> to the best of my ability to measure distances on these ~archD-sized
> maps).  Looking at the 10000 meter grid-lines on the quad, they are indeed
> not parallel to the edges of the plotted box (nor to the edges of the paper).
> These observations suggest to me that the 100k quad maps are not actually
> plotted with a UTM projection. However, when reading the other map information,
> I cannot see other projection information, other than (perhaps) something to do
> with aerial photograph planimetry.
>
> The only two -J projections for geographic data that I know off the top of my head that would
> make lon/lat lines be perfectly parallel to the edges of a map are: -JM|m and -JX|x...d
> The equatorial Mercator option, -Jm1:100000  gives a decent overlay, but it seems funny to me
> that the topo-map producers would use a projection with such distortion far away from the
> Equator.  If my target plot had any more latitude range to it, I suspect that -Jm would overlay
> unsuitably.
>
> For the fake-it-as-a-linear-plot-with-the-R-coordinates-in-degrees option, -Jx1:100000d, I get
> the -J syntax error message:
> % psbasemap -R-120:00/-119:00/34:30/35:00 -Jx1:100000d -Ba0.5 -X1 -Y1 --PAPER_MEDIA=archD >! test.Jx.ps
> psbasemap: GMT SYNTAX ERROR -J option.  Correct syntax:
> ....rest of the error message...
>
> Two questions:
> 1) Has anyone perfectly-overlain GMT-produced *.ps files onto 1:100k
> quads?  If so, with what -J projection option.  Some searching through
> this list and the whole Internet didn't produce wisdom on this.  I looked
> at a few 1:24,000 USGS topo quads that I have, and these seem to
> have the same issue: the map legend indicates a polyconic projection
> but the boundary of the maps are perfectly at the 7.5-minute lon-lat
> values and these map-boundaries are perfectly parallel to the edges
> of the printed map.
>
> 2) Should -Jx1:100000d produce a usage error?  -JX1d does not produce
> an error, and it yields a plot with FANCY basemap type as I would
> expect for a geographic-aware map.
>
>
> By the way, the GeoTIFF associated with this quad has metadata (tags) indicating
> that this file is indeed in UTM 11 North; when I view this in OpenEV, the lon-lat
> edges of the map's outer box are indeed not parallel to the graphic window,
> and the skewed-angle looks just like that of the psbasemap -Ju+11... map.
>
> I may ask my ESRI-savvy colleagues how they would produce
> an overlay map, but thought I'd try this list first.
>
> Thanks!
>
>
>
>
> To unsubscribe, send the message "signoff gmt-help" to listserv@...
>

To unsubscribe, send the message "signoff gmt-help" to listserv@...

[Paul Wessel-2]

Dan-

There is a few bugs in -Jx1:<scale>.  I am working on it.
-p

On Oct 6, 2009, at 7:54 AM, Dan Scheirer wrote:

> Greetings,
>
> I would like, but have not been able, to plot data in GMT so that I  
> can superimpose
> my PS file onto digital (PDF) versions of USGS 1:100,000 (60'x30')
> topographic quadrangle maps.  The quad in question is in central
> California (Cuyama), at the western edge of UTM zone 11 North,
> but I think that my question would apply to many/all such
> 100k quad maps produced by the USGS (and perhaps the other
> scale series, too).
>
> If you have experience/interest with this, please read-on...
>
> Thanks,
> Dan Scheirer
>
>
> I thought this would be easy to do because the map
> legend says: "Projection and 10000 meter grid, zone 11
> Universal Transverse Mercator".  I tried:
> psbasemap -Ju+11/1:100000 -Ba0.5 -R-120/-119/34.5/35 -X1i -Y1i >!  
> test.ps
> It's immediately apparent that the lon-lat grid-lines are not parallel
> to the edges of the paper, esp. because this is on the edge of a
> UTM zone, and this makes sense to me.  If I use the -R<lower-left>/
> <upper-right>r
> specification, I get edges to the map-box parallel to the edges of  
> the paper, but
> I do not get the full lat-lon range as specified above (thin  
> "triangles" are missing).
> This also makes sense to me given the TM distortion, maximal near  
> the edge of a zone.
>
> What is unexpected to me is that the USGS quad map has the min-max  
> longitudes
> and min-max latitudes perfectly define the corners of the map, and  
> that
> these corners define edges that are parallel to the edges of the  
> paper (parallel
> to the best of my ability to measure distances on these ~archD-sized
> maps).  Looking at the 10000 meter grid-lines on the quad, they are  
> indeed
> not parallel to the edges of the plotted box (nor to the edges of  
> the paper).
> These observations suggest to me that the 100k quad maps are not  
> actually
> plotted with a UTM projection. However, when reading the other map  
> information,
> I cannot see other projection information, other than (perhaps)  
> something to do
> with aerial photograph planimetry.
>
> The only two -J projections for geographic data that I know off the  
> top of my head that would
> make lon/lat lines be perfectly parallel to the edges of a map are: -
> JM|m and -JX|x...d
> The equatorial Mercator option, -Jm1:100000  gives a decent overlay,  
> but it seems funny to me
> that the topo-map producers would use a projection with such  
> distortion far away from the
> Equator.  If my target plot had any more latitude range to it, I  
> suspect that -Jm would overlay
> unsuitably.
>
> For the fake-it-as-a-linear-plot-with-the-R-coordinates-in-degrees  
> option, -Jx1:100000d, I get
> the -J syntax error message:
> % psbasemap -R-120:00/-119:00/34:30/35:00 -Jx1:100000d -Ba0.5 -X1 -
> Y1 --PAPER_MEDIA=archD >! test.Jx.ps
> psbasemap: GMT SYNTAX ERROR -J option.  Correct syntax:
> ....rest of the error message...
>
> Two questions:
> 1) Has anyone perfectly-overlain GMT-produced *.ps files onto 1:100k
> quads?  If so, with what -J projection option.  Some searching through
> this list and the whole Internet didn't produce wisdom on this.  I  
> looked
> at a few 1:24,000 USGS topo quads that I have, and these seem to
> have the same issue: the map legend indicates a polyconic projection
> but the boundary of the maps are perfectly at the 7.5-minute lon-lat
> values and these map-boundaries are perfectly parallel to the edges
> of the printed map.
>
> 2) Should -Jx1:100000d produce a usage error?  -JX1d does not produce
> an error, and it yields a plot with FANCY basemap type as I would
> expect for a geographic-aware map.
>
>
> By the way, the GeoTIFF associated with this quad has metadata  
> (tags) indicating
> that this file is indeed in UTM 11 North; when I view this in  
> OpenEV, the lon-lat
> edges of the map's outer box are indeed not parallel to the graphic  
> window,
> and the skewed-angle looks just like that of the psbasemap -Ju+11...  
> map.
>
> I may ask my ESRI-savvy colleagues how they would produce
> an overlay map, but thought I'd try this list first.
>
> Thanks!
>
>
>
>
> To unsubscribe, send the message "signoff gmt-help" to listserv@...

To unsubscribe, send the message "signoff gmt-help" to listserv@...

[Dan Scheirer]

Greetings,

A resolution to my initial question of a few weeks ago, plus a mini-bug report.

Thanks,
Dan



It turns out that the map-makers of the 1:100k quadrangle in question used a transverse
Mercator projection with a local origin (at the bottom, middle of the quadrangle)
when they constructed the "paper" map that was downloadable in the PDF file.
The parameters of that local origin were viewable in the *meta.txt metadata
file associated with this publication.  Curiously, they were embedded in a
section entitled "UTM Zone 11", and I did not see them initially.  All of the
digital coordinates of the GIS files were indeed in UTM Zone 11 meters, so the local origin
was only used to make the "paper" map less skewed.  It also turns out that if I
expanded the PDF file greatly, I could see the local transverse Mercator skew in
slightly-diverging east-west boundaries of the quadrangle and slightly concave-down
north-south boundaries.  No information about the local transverse Mercator
projection was imprinted on the map (seems like a significant omission to me);
it only referred to UTM Zone 11.

The plus-side of making a quadrangle with individually-selected local projections
is that the lon-lat lines are very-nearly-parallel to the edges of the paper.  The down-side
of making a quadrangles in this way is that neat-lines between adjacent maps
won't match perfectly.  In talking with a number of USGS colleague map-makers
(I work at the USGS), I gather that there are several ways to make "paper" maps
less-skewed than would be implied by using the projected coordinate system
intrinsic to the underlying digital data.  I'm not sure if people care about neat-lines
any more...

Mini-bug, discovered while trying the map-matching as described above:

I believe that the -R region-checking could be improved for cases with
non-geographic projections (e.g. -JX).  Specifically, when using
-R<xmin>/<ymin>/<xmax>/<ymax>r, there is not test that xmin<xmax
and ymin<ymax.  This test is present for all projections that do not
use the -R...r specification.  The following 5 commands illustrate this behavior;
the first 2 make maps as expected; the latter 3 do not generate
errors even though the -R values are inconsistent, and they yield
maps that do not make sense.

psbasemap -JX5i -Ba1 -R0/5/21/26 -P -X1.5i -Y2i >! test.JX.no_r.ps
psbasemap -JX5i -Ba1 -R0/21/5/26r -P -X1.5i -Y2i >! test.JX.good_r.ps
psbasemap -JX5i -Ba1 -R0/26/5/21r -P -X1.5i -Y2i >! test.JX.bad1_r.ps
psbasemap -JX5i -Ba1 -R5/26/0/21r -P -X1.5i -Y2i >! test.JX.bad2_r.ps
psbasemap -JX5i -Ba1 -R5/21/0/26r -P -X1.5i -Y2i >! test.JX.bad3_r.ps

I'm not sure whether or not the errors above should also apply to geographic
projections.  The following 5 lines illustrate very similar behavior with -JM:
no error messages and non-sensical maps.

psbasemap -JM5i -Ba1 -R0/5/21/26 -P -X1.5i -Y2i >! test.JM.no_r.ps
psbasemap -JM5i -Ba1 -R0/21/5/26r -P -X1.5i -Y2i >! test.JM.good_r.ps
psbasemap -JM5i -Ba1 -R0/26/5/21r -P -X1.5i -Y2i >! test.JM.bad1_r.ps
psbasemap -JM5i -Ba1 -R5/26/0/21r -P -X1.5i -Y2i >! test.JM.bad2_r.ps
psbasemap -JM5i -Ba1 -R5/21/0/26r -P -X1.5i -Y2i >! test.JM.bad3_r.ps

For the "incorrect" -R...r specifications, it's possible that they can
correctly refer to most of the rest-of-the-world-not-included-in-the-
non-R...r-specification, and perhaps the bad-maps arise simply
because -JM cannot deal with this.



----- Original Message ----

> From: Mike Willis <willismi@...>
> To: gmt-help@...
> Sent: Mon, October 5, 2009 2:00:47 PM
> Subject: Re: [GMT-HELP] how to overlay psbasemap onto USGS 1:100k quads?
>
> If you scanned the map, then the image is already projected, so you
> want the linear scale. Take a look at example 28.
>
> Cheers
>
> On Mon, Oct 5, 2009 at 4:54 PM, Dan Scheirer wrote:
> > Greetings,
> >
> > I would like, but have not been able, to plot data in GMT so that I can
> superimpose
> > my PS file onto digital (PDF) versions of USGS 1:100,000 (60'x30')
> > topographic quadrangle maps.  The quad in question is in central
> > California (Cuyama), at the western edge of UTM zone 11 North,
> > but I think that my question would apply to many/all such
> > 100k quad maps produced by the USGS (and perhaps the other
> > scale series, too).
> >
> > If you have experience/interest with this, please read-on...
> >
> > Thanks,
> > Dan Scheirer
> >
> >
> > I thought this would be easy to do because the map
> > legend says: "Projection and 10000 meter grid, zone 11
> > Universal Transverse Mercator".  I tried:
> > psbasemap -Ju+11/1:100000 -Ba0.5 -R-120/-119/34.5/35 -X1i -Y1i >! test.ps
> > It's immediately apparent that the lon-lat grid-lines are not parallel
> > to the edges of the paper, esp. because this is on the edge of a
> > UTM zone, and this makes sense to me.  If I use the
> -R/r
> > specification, I get edges to the map-box parallel to the edges of the paper,
> but
> > I do not get the full lat-lon range as specified above (thin "triangles" are
> missing).
> > This also makes sense to me given the TM distortion, maximal near the edge of
> a zone.
> >
> > What is unexpected to me is that the USGS quad map has the min-max longitudes
> > and min-max latitudes perfectly define the corners of the map, and that
> > these corners define edges that are parallel to the edges of the paper
> (parallel
> > to the best of my ability to measure distances on these ~archD-sized
> > maps).  Looking at the 10000 meter grid-lines on the quad, they are indeed
> > not parallel to the edges of the plotted box (nor to the edges of the paper).
> > These observations suggest to me that the 100k quad maps are not actually
> > plotted with a UTM projection. However, when reading the other map
> information,
> > I cannot see other projection information, other than (perhaps) something to
> do
> > with aerial photograph planimetry.
> >
> > The only two -J projections for geographic data that I know off the top of my
> head that would
> > make lon/lat lines be perfectly parallel to the edges of a map are: -JM|m and
> -JX|x...d
> > The equatorial Mercator option, -Jm1:100000  gives a decent overlay, but it
> seems funny to me
> > that the topo-map producers would use a projection with such distortion far
> away from the
> > Equator.  If my target plot had any more latitude range to it, I suspect that
> -Jm would overlay
> > unsuitably.
> >
> > For the fake-it-as-a-linear-plot-with-the-R-coordinates-in-degrees option,
> -Jx1:100000d, I get
> > the -J syntax error message:
> > % psbasemap -R-120:00/-119:00/34:30/35:00 -Jx1:100000d -Ba0.5 -X1 -Y1
> --PAPER_MEDIA=archD >! test.Jx.ps
> > psbasemap: GMT SYNTAX ERROR -J option.  Correct syntax:
> > ....rest of the error message...
> >
> > Two questions:
> > 1) Has anyone perfectly-overlain GMT-produced *.ps files onto 1:100k
> > quads?  If so, with what -J projection option.  Some searching through
> > this list and the whole Internet didn't produce wisdom on this.  I looked
> > at a few 1:24,000 USGS topo quads that I have, and these seem to
> > have the same issue: the map legend indicates a polyconic projection
> > but the boundary of the maps are perfectly at the 7.5-minute lon-lat
> > values and these map-boundaries are perfectly parallel to the edges
> > of the printed map.
> >
> > 2) Should -Jx1:100000d produce a usage error?  -JX1d does not produce
> > an error, and it yields a plot with FANCY basemap type as I would
> > expect for a geographic-aware map.
> >
> >
> > By the way, the GeoTIFF associated with this quad has metadata (tags)
> indicating
> > that this file is indeed in UTM 11 North; when I view this in OpenEV, the
> lon-lat
> > edges of the map's outer box are indeed not parallel to the graphic window,
> > and the skewed-angle looks just like that of the psbasemap -Ju+11... map.
> >
> > I may ask my ESRI-savvy colleagues how they would produce
> > an overlay map, but thought I'd try this list first.
> >
> > Thanks!
> >
> >
> >
> >
> > To unsubscribe, send the message "signoff gmt-help" to listserv@...
> >
>
> To unsubscribe, send the message "signoff gmt-help" to listserv@...

     

To unsubscribe, send the message "signoff gmt-help" to listserv@...