Historical USGS Quads

This will be a very quick post today.  We wanted to use some historical USGS quads for analysis.  We had a ready source, used ArcGIS to georeference them to a grid, but then wanted something approaching a seamless product.

Enter gdal.

First step, create an image to put the mosaic into:


gdal_merge -o usgs_merge.tif -createonly input1.tif input2.tif ...


Then mosaic into that.  Since we have a shapefile with the shape of the bounds of each of the images in question, we can use a cutline to perform the mosaic, e.g.

gdalwarp -cutline quadrangle_index_24k.shp -cwhere "quadname_ = West_View" "West_View.tif" usgs_merge.tif

In all honesty, this was scripted in Excel for expediency, but could easily be done in Windows Command with for-in-do, or BASH at a UNIX machine.

Final result:

Geometry Collections and Small Headaches (pics)

Geometry Collections and Small Headaches

The disk that held the personal geodatabases of our contour datasets died a while back, but not before I loaded the contours into PostGIS and started serving them up. Our new intern is working on putting together some shapefiles layer groups in ArcGIS for map production, and asked for a missing one… . I don’t have it, but in principle can always extract it from PostGIS and serve it back as a shapefile.

So I performed an intersection:

CREATE TABLE public.contour_2_cut AS
	SELECT contour.elevation, (ST_Intersection(contour.the_geom), cutline) AS the_geom, gid FROM
		public.contour_2 AS contour,  (ST_GeomFromText('POLYGON((2000000 860000, 2500000 860000, 2500000 460000, 2000000 460000,2000000 860000))', 3734)) AS cutline;

After adding a primary key:

ALTER TABLE public.contour_2_cut ADD PRIMARY KEY (gid);

I couldn’t display it in uDig, which complained that it was a geometry collection. Ah, of course. In the ESRI world (whence I hark) an intersection takes the form of the least of the geometries involved, so the intersection of a polyline and a polygon is a polyline. Not so in PostGIS. PostGIS supports geometry collections, so naturally the intersection includes both types of geometry. (Side note– the viewing in uDig didn’t fail until it got to the polygon– it first loaded all the contour polylines, then discovering the polygon, threw an error).

What to do? Well, I hoped pgsql2shp would bail me out:

pgsql2shp -f contours_block_08_01 -h cmac-srv-gis -u postgres CM public.contour_2_cut
Initializing... type 'GEOMETRYCOLLECTION' is not Supported at this time.
The DBF file will be created but not the shx or shp files.
You've found a bug! (pgsql2shp.c:2864)

But it doesn’t support GEOMETRYCOLLECTIONs either.

Time to dig deeper. So how do we tell what a geometry is?

SELECT ST_GeometryType(st_geom) FROM contour_2_cut;

and I see ST_LineString is the descriptor for our lines. So, let’s create a copy of the table with just linestrings (we could remove records from the existing table as well):

CREATE TABLE contour_2_cut_line AS
	SELECT * FROM contour_2_cut WHERE
		ST_GeometryType(the_geom) ='ST_LineString';

Parcel Annotations in GeoServer (with some Maplex help) (cont. 2)

I promised pics from our labeled parcels:

 

http://smathermather.wordpress.com/2011/02/01/parcel-annotations-in-geoserver-with-some-maplex-help/

http://smathermather.wordpress.com/2011/02/01/parcel-annotations-in-geoserver-with-some-maplex-help-cont-1/

Parcel Annotations in GeoServer (with some Maplex help)

smathermather:

We have a guest blogger today– Ramon, a bright and hard-working intern I’ve had the pleasure of working with for over a year.  If you’re looking for someone versed in Postgre/PostGIS/GeoServer/OpenLayers, Ramon’s been my rock as we’ve been building our internal system, doing everything from basic grunt work to esoteric trouble-shooting.  I’m trying hard not to fall on my face now that he’s moved on.

We wanted to take advantage of the advance labeling of Maplex for parcels in GeoServer.  GeoServer SLDs don’t yet (I think) have automatic rotation for best fit for polygons, but Maplex (an ESRI label extension that’s rolled in with the ArcINFO license in ArcGIS desktop) does.  Here’s how we took advantage of that:

Ramon:
The general procedure to generate the attribute table is as follows: (Note: It is best to do this in sections because in the limitations of labels that could be generated per annotation feature)
1) Generate the labels using Maplex as the label engine in Arcmap.
2) Convert the labels into annotation (in a database)
3) Join the annotation table back to the original polygon that needs labeled.
4) Export the dataset to create a new shapefile with the combined attribute table of the original shapefile and “annotation feature”
5) Import the shapefile to PostgreSQL

Notes: I added the following fields to the shapefie before importing it to PostgreSQL
a) acreage – area of the parcel in acres
b) res_prop – “Y” if parcel is within 1500 ft of our area of interest, otherwise “N”
c) ang – small integer conversion of the “Angle” value generated by Maplex (this may be scrapped from the workflow in the future)
d) Rot_Ang – the rotation angle in terms of Geoserver convention.
- calculated as “zero” minus “ang” (see the UPDATE SQL below)
e.g. UPDATE base.parcel_annotations_med SET rot_ang = 0 – ang;

smathermather:

This final step step converts annotation rotation values, which are rotated from horizontal up to 90 degrees either in a positive or negative to GeoServer label rotation convention, which run the full arc of a circle.

We’ll have a post that follows with the final SLD, and screen shot of the great labeling effect.  Stay posted.

Garmin Custom Maps

I don’t spend much time on field equipment, e.g. GPS units.  Why?  Lack of interest?  No.  Lack of space in my job description.  Yes, I know.  A GIS Manager who doesn’t have a GPS is it’s own category of silly.  None-the-less, when a colleague told me that Garmin now has firmware support for custom maps for some of it’s GPS units, I told him it might be a good choice for his field work and that I’d work on loading our maps onto his unit when it came in.

Sure enough, it came in, and we’ve loaded them.  You can follow Garmin’s directions, but if you already have georeferenced data, skip them and by whatever means you know, generate a jpeg with a world file (for you ArcMap users, goto data view, file:export choose jpeg, and check the world file checkbox).  You can export an image from any desktop or server package, and if you know the upper left coordinate of the image, and the pixel size, you have enough to create a world file.  Easiest hand crafting of geographic data you can imagine (it’s the gateway for bigger and better hand crafted geographic data).

Now, you can use MapTiler to generate the kind of KMZ that is so near and dear to the garmin’s heart.  Just choose KMZ as your output type, choose one of the garmin types in the next screen, and the rest easy easy GUI fun.

Limitations:

Well, Garmin doesn’t want you to load more than 100 tiles (probably for performance reasons), but what does that mean practically?  For the testing I did, that meant for a 4 foot pixel (most of our best data isn’t more spatially accurate than that, even if it’s more precise), I was able to load about 3200 acres– enough for one of our larger parks.  Good enough for our work, but you might consider a deeper tool if you need more raster data than that at a time.

ArcGIS, Layer Packages and Workarounds

I’ve said before that I work in a hybrid shop– part ESRI, part Open Source. We use GeoExt/OpenLayers/GeoServer/PostGIS for the Enterprise stuff, and ESRI for the analyses, some specialized cartography, and the easy way out on one-off projects. I mostly blog about the Open Source stuff, ’cause that’s where my heart (and budget) lays.

So, now for some ESRI Fan geek-dum. Plus the obligatory work around. Fair warning:

Probably my favorite feature introduced with 9.3 is Layer Packages. Right click on a layer group, choose to export to layer package, and ArcGIS packages up your layer as a compressed layer group that you can send to friends, colleagues, and enemies.

Problem– sometimes the layer groups don’t work when you try to open them at the other end. So I’ve taken to testing them– create, double click to test (and thus decompress) and make sure there are no memory errors part-way through. I discovered a good work-around when there are memory errors– the packages decompress just fine with 7-Zip. What I’ve taken to doing now to ensure my friends, colleagues, and enemies can use my layer packages is to create, unzip with 7-zip, and re-zip with 7-zip as a *.zip file. 7-Zip and WinZip, etc. are much more stable than the engine used by Esri, and I haven’t had errors in the creation of the layer group, just the automatic decompression of it.

In essence, I’m using layer packages to help me organize the data to send, but using good old regular compression software to deliver it.

Shapfile spatial indices and PostGIS dumps

Another short post… . I have in place a framework for maintenance of our boundary data that works really well. I feed an edited shapefile into the database, and PostGIS creates all the versions of the boundary that need created. The database then exports the managed data to shapefile so there is a database version and a shapefile version depending on need.

When all was said and done, I had several shapefiles that displayed at some zoom levels in ArcGIS, but not others. Google to the rescue:

http://forums.esri.com/Thread.asp?c=93&f=1730&t=164942

I had outdated spatial indices that are not created by pgsql2shp, and since they didn’t reflect the real geometry, ArcGIS was choking on them and not displaying them properly. In the above link, Kirk Kuykendall suggests fixing the index in ArcCatalog. For me, the geometry is simple enough, I just deleted the index (*.sbn) with no noticeable performance hit.