Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Posts Tagged ‘ArcGIS’

Proper (ab)use of a database, contour interpolation using #postgresql #postgis #arcgis

Posted by smathermather on August 6, 2012

Anyone who has been following along at home knows I don’t think much like a DBA.  Sometimes that’s good; mostly it’s probably bad.  In this post, I hope it will be interesting.

The problem of the day is how to take engineering contours derived from breaklines, a lidar point cloud, and all the lot, and do a good job interpolating that to a DEM.  This is an interesting problem, as it’s another one of those problems with so many sub-optimal solutions.  In part, it’s an impossible problem– going from a surface model to contours is inherently a loss of information, so going the other way is difficult to do believably, both with respect to believable smoothness and the requisite feature retention.  This is a little like the pugilistic pie puppets problem: ‘flaky vs. tender’ of making good pie crust.  Doing one well seems to preclude doing the other well.

Now, Chris Hormann, of PovRay fame has a solution, there are solutions in ArcGIS, and in GRASS, and plenty of other attempts.  Here’s a decent list, which I haven’t checked for timeliness:  http://vterrain.org/Elevation/contour.html.  Since we have a ArcGIS license, it handles the dataset I have OK, the results checked out OK, this is what we chose to use for now.  I tried GRASS, but it was too slow for the amount of data being processed.  The command in ArcGIS ends up looking something like this:


TopoToRaster_sa "N2260635.shp elevation Contour;N2260640.shp elevation Contour;N2260645.shp elevation Contour;N2265633.shp elevation Contour;N2265640.shp elevation Contour;N2265645.shp elevation Contour;N2270635.shp elevation Contour;N2270640.shp elevation Contour;N2270645.shp elevation Contour;N2275650.shp elevation Contour" "C:\Temp\TopoToR_N2269.tif" 2 "2265031 645070 2270065 639945" 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # # #

And thus we go from this:

To this:

Now to do it 623 more times… .  And in walks PostgreSQL for some abuse as a spatial engine with string manipulation to boot.

Each tile of contours I run through this approach needs the adjacent 8 tiles of contours loaded in order to avoid dreaded edge effects:

Except that sometimes, there aren’t 8 tiles available, like on the edge of large water bodies:

So, we can’t just use a heuristic that loads the 8 adjacent tiles– we need to enumerate these tiles.  An ST_Touches should do the trick here:

SELECT p.tile, a.tile AS tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom);

But, that does not arrange our table very nicely:

I want something that looks like this:

N2225655, N2225660, N2220650 …

N2225660, N2225665, N220665 … etc.

Thank goodness for Postgres’ String_agg, and moreover for Leo Hsu and Regina Obe’s blog.

Adapting their entry, we can reorganize those entries like this:

SELECT
STRING_AGG(a.tile, ',' ORDER BY a.tile)
As tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom)
GROUP BY p.tile
ORDER BY p.tile;

Ah, now that’s better:

And if we want, we can use an arbitrary delimiter, and build our entire command in one swell foop:

SELECT p.tile,
'TopoToRaster_sa "'
|| 'W:\contours\cuy_contours\shapefiles\' || p.tile || '.shp elevation Contour;'
|| 'W:\contours\cuy_contours\shapefiles\'
|| STRING_AGG(a.tile, '.shp elevation Contour;W:\contours\cuy_contours\shapefiles\' ORDER BY a.tile)
|| '.shp elevation Contour;" '
|| '"v:\svm\dem\cuyahoga\' || p.tile || '.tif" 2'
|| ' "' || replace(replace(replace(ST_Extent(p.geom)::text, 'BOX(', ''), ')', ''), ',', ' ') || '"'
|| ' 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # '
As tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom)
GROUP BY p.tile
ORDER BY p.tile

Ok.  I have 623 more tiles to do.  Until next time… .

Posted in Analysis, ArcGIS, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , , , , , , , , | Leave a Comment »

Historical USGS Quads

Posted by smathermather on March 4, 2011

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:

Posted in Other | Tagged: , , , | 1 Comment »

Geometry Collections and Small Headaches (pics)

Posted by smathermather on February 10, 2011

Posted in Database, GeoServer, PostGIS, SQL | Tagged: , , , , | Leave a Comment »

Geometry Collections and Small Headaches

Posted by smathermather on February 10, 2011

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';

Posted in Database, PostGIS, SQL | Tagged: , , , , | 5 Comments »

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

Posted by smathermather on February 4, 2011

I promised pics from our labeled parcels:

 

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

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

Posted in GeoServer, PostGIS | Tagged: , , , , , , , , | Leave a Comment »

Parcel Annotations in GeoServer (with some Maplex help)

Posted by smathermather on February 1, 2011

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.

Posted in GeoServer, PostGIS | Tagged: , , , , , , , , | 2 Comments »

Garmin Custom Maps

Posted by smathermather on January 11, 2011

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.

Posted in Uncategorized | Tagged: , , , , | 3 Comments »

ArcGIS, Layer Packages and Workarounds

Posted by smathermather on December 22, 2010

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.

Posted in Other | Tagged: | 2 Comments »

Shapfile spatial indices and PostGIS dumps

Posted by smathermather on July 19, 2010

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.

Posted in Database, PostGIS, SQL | Tagged: , , , | 2 Comments »