Kurt's Weblog

vesseltrax is a nice website with AIS ship tracking of the Galveston/Houston port area with Google Maps and Earth based interfaces. The google earth file has nicer ship icons than I have.

The front page has a drop down that lets you get the Google Earth kml file.



Here is the Google Maps interface that shows a ships AIS message 5 name and cargo information.



The Google Earth view with little ship icons.

I am still trying to decide on what database packages if any a fink sqlobject package should depend, but I have a draft sqlobject-py.info file that you can put in /sw/fink/10.4/local/main/finkinfo/databases

Note, that if you are readying this long after I wrote this post, the sqlobject will hopefully be in the fink tree and not on my drafts page.

My draft fink packages - download sqlobject-py.info dated 31-Jan-2007 or newer.

BTW, if someone really knows sqlalchemy and/or sqlobject, I could really use a tutorial on how to use these database wrappers effectively.

I've been trying to recapture that missed shot. No luck on the lighting this morning, but I did capture this shot of the forest.

I don't have any details since this is a fairplay headline:

Navy sub commander relieved of duty
.
THE commander of a US nuclear submarine has been relieved of duty
following his boat's collision with a tanker in the Persian Gulf.

Video Tutorials for SketchUp to Google Earth [googleearthblog.com]

Shape2Earth - GIS Data to Google Earth talks about converting shape files to Google Earth KML using some commercial windows only tools. Here is another option. I've shown this before with s57 ENC chart data, but here it is with a Arc shape file. (Using gdal version 1.4.0 or newer)

 ogrinfo sbnms.shp
INFO: Open of `sbnms.shp'
      using driver `ESRI Shapefile' successful.
1: sbnms (Polygon)

So there should be a polygon in that file.

ogr2ogr -f KML sbnms.kml sbnms.shp

Which creates this kml

<?xml version="1.0" encoding="utf-8" ?>
<kml xmlns="http://earth.google.com/kml/2.0">
<Document><Folder><name>sbnms</name>
  <Placemark>
    <description><![CDATA[
      <b>OBJECTID:</b> <i>1</i><br />
      <b>AREA:</b> <i>0.23956600000</i><br />
      <b>PERIMETER:</b> <i>2.04266000000</i><br />
      <b>Shape_Leng:</b> <i>2.04266456234</i><br />
      <b>Shape_Area:</b> <i>0.23956563571</i><br />
   ]]></description>      <Polygon><outerBoundaryIs><LinearRing><coordinates>-70.502624517077535,42.651229858376865 -70.216651917737892,42.766723627364335 -70.035064690795366,42.093299872372249 -70.054344171734087,42.1051712037602 -70.067070010211893,42.108303073944938 -70.086586000518764,42.117519381620241 -70.106071466728224,42.120601658877817 -70.123886112296802,42.126754758125287 -70.140045164348635,42.128528592356282 -70.154968262866575,42.133419041708223 -70.172912604417078,42.13481140155551 -70.196044924942854,42.132099157383799 -70.217071529781094,42.133392331376747 -70.238891607484177,42.129699711869044 -70.2558517459629,42.124347685078639 -70.277999877386733,42.115261073385255 -70.470436094985956,42.12923813443058 -70.597366328254395,42.551074986184581 -70.586975101270227,42.558498376829306 -70.58387756894885,42.563468929236009 -70.572532649020488,42.575222019909717 -70.555580138073651,42.580753321639378 -70.541786193673204,42.587898251490522 -70.528427130349584,42.595039367575716 -70.51586913948303,42.606506341984264 -70.50587462797651,42.621070862684029 -70.501319887816891,42.633117669215117 -70.501304632753104,42.642444604042602 -70.502624517077535,42.651229858376865</coordinates></LinearRing></outerBoundaryIs></Polygon>
  <Style><LineStyle><color>ff0000ff</color></LineStyle>  <PolyStyle><fill>0</fill></PolyStyle></Style>
  </Placemark>
</Folder></Document></kml>

The results in Google Earth look like this:

Princess Cruise Lines pays fine, restitution in whale's death

by Jill Burke:
Princess Cruise Lines is facing a significant fine after one of its
ships struck and killed a humpback whale in Southeast Alaska.



The collision took place in July 2001 near Glacier Bay National Park
and Preserve.



Princess Cruise Lines is expected to enter a guilty plea Monday in
U.S. district court in Anchorage. The company is also expected to pay
a $755,000 fine because the ship Dawn Princess failed to operate at a
slow speed around the whale.

SISTER SANCTUARIES TO PROTECT ENDANGERED WHALES AT BOTH ENDS OF ANNUAL MIGRATION [noaa.gov]

United States and Dominican Republic Partner in Historic Conservation Effort

NOAA image of humpback whale breaching in the NOAA Stellwagen Bank
National Marine Sanctuary, which is one of the acrobatic behaviors
that whale watchers enjoy seeing.Jan. 29, 2007 - NOAA established a
"sister sanctuary" arrangement between the NOAA Stellwagen Bank
National Marine Sanctuary off the coast of Massachusetts and the
Marine Mammal Sanctuary of the Dominican Republic, two marine
protected areas 3,000 miles apart that provide conservation programs
for the same population of humpback whales.

Several of us were discussing information overload for mariners in the last day. It is always tough to know how and how much information to display to users.

T.M.I., Robo-Dude [defensetech.org]

But that's too many robots, Hudas says. Four drones is the realistic
max. And a third crewman at an additional console is ideal. And that's
assuming a minimal level of human intervention in the drones
activities. Basically, you tell a drone what to do, confirm the
command, then let it go. Now, if the drone wants to kill something,
it's going to need a soldier's permission. But for surveillance and
reconnaissance, it can make its own decisions. "With those
applications," Hudas says, "we don't even want a soldier."

Find out who is where. Hard to dig up ships unless you know what ports they are using, but interesting. I don't see a how it works on their web page. I would be interested in knowning which databases they are scraping.

Just a reminder that the NSF visualization competion is open for submissions until May 31st.

NSF 2007 Science and Engineering Visualization Challenge

There was a light snow last night. Here is what Gregg Hall looks like from the woods.

I am working my way towards having PCL be in fink for mac users. I just checked in quadtree-py and owslib-py which are both useful in their own rights.

Quadtree

OWSLib provides OWS/WMS/WMF capabilities (OCG Web Services)

I need to export all the frames in about 10 HD sized AVI files so that I can overlay figures and text. I tried Quicktime Pro's export functionality. File - Export - Movie to Image Sequence - Options - PNG. That took about an hour and a half to export a 3:20 minute avi. Trying mplayer with altivec on fink looks to take under 30 minutes. Much better.

mplayer -nosound -vo png z_gray.avi

activecaptain.com

Requires login to use. It is a Yahoo Maps based application. Unfortunately, if you submit information to this site, you can only use that info from within the site.

Yesterday, we had a talk by Suzanne Carbotte of LDEO about access to marine datasets.

Serving Multidisciplinary Data for Marine Geoscience Research:
     Approaches and Lessons Learned
Abstract: Digital data collections serve an increasingly important
role in research and education: they aid data discovery, democratize
access to research resources, enable verification of research results,
and provide new research opportunities. The Marine Geoscience Data
System, hosted at LDEO, is a resource collection that serves
multi-disciplinary marine geoscience data for the NSF Ridge2000 and
MARGINS programs and for Southern Ocean Research conducted on the RV
Palmer.



An overview of system architecture including data and access services
will be presented. Our experience with this system highlights issues
associated with managing such collections, including the need to
embrace rapidly evolving technologies, the need for new approaches to
metadata collection and data transfer to archives, and the importance
of close interaction with user

I am excited that she put so much emphasis on WMS/WMF. She also talked about OAI-PMH. Not stuff I had heard of before. That stands for "Open Archives Initiative - Protocol for Metadata Harvesting." openarchives.org/

Standards for Repository Interoperability
The Open Archives Initiative develops and promotes interoperability
standards that aim to facilitate the efficient dissemination of
content. OAI has its roots in the open access and institutional
repository movements. Continued support of this work remains a
cornerstone of the Open Archives program. Over time, however, the work
of OAI has expanded to promote broad access to digital resources for
eScholarship, eLearning, and eScience.

This morning's sunrise had people running for their cameras - me included. I just missed the steam coming out of the stack looking like it was on fire. Just missed it by less than a minute. Here is what I did capture, which is not bad, but timing is everything.

This is one of those basic problems. How do you know if a point is within a polygon. Yes, I have coded this up myself, but I would rather use a solid package to do this and get lots of other operators done right. I just got PCL-Core working and have a little example that shows testing points against a polygon.

Note that WKT is the "Well Known Text" form of geographic data. I seem to be unable to get it right when typing. I think my history with World ToolKit (WTK) is messing with my head.

Here is the example:

In [1]: from cartography.geometry import Geometry, Point, LinearRing, Polygon
In [2]: wkt = 'POLYGON ((-1.0 50.5, -0.5 51.2, 0.3 50.9, -1 50.5))'
In [3]: a = Geometry.fromWKT(wkt)
In [4]: wtk_POINT_IN = 'POINT (-0.5 51.15)'
In [5]: pIn = Geometry.fromWKT(wtk_POINT_IN)
In [6]: wtk_POINT_IN = 'POINT (-0.5 51.5)'
.
In [7]: pIn.within(a)
Out[7]: 1
.
In [8]: wtk_POINT_OUT = 'POINT (100 51.5)'
In [9]: pOut = Geometry.fromWKT(wtk_POINT_OUT)
.
In [10]: pOut.within(a)
Out[10]: 0

That's it. I am still working on the fink packaging (which is giving me a little bit of trouble), but should have packages up soon. http://vislab-ccom.unh.edu/~schwehr/software/fink/. You will need owslib-py.info, quadtree-py.info and a couple pcl-*.info files.

The temperature today is a nice warm -2 F. I captured this image on a window today: Ice crystals growing on a window pane.

periodic table of visualization techniques

Long Beach energy project halted [latimes.com]

...
In a unanimous vote, the Long Beach Board of Harbor Commissioners
decided to end an environmental review of the project that was
launched more than two years ago but had slipped far behind
schedule. The action effectively terminates the effort by the port and
a partnership of Mitsubishi Corp. and ConocoPhillips to build a
$700-million liquefied natural gas plant inside the busiest cargo port
in the nation.



"The project is dead," said Doris Topsy-Elvord, a commissioner and
former Long Beach councilwoman. "It's been a long haul, we've worked
hard at it, but it is at a standstill. I do not think there's a
possibility it will come back," she said.
...
Since Long Beach officials began reviewing the project, other energy
companies have proposed liquid natural gas terminals offshore in
Southern California where the risk to the public is minimal.



Australia-based BHP Billiton proposes a floating terminal off the
coast of Ventura County; Crystal Energy Co. wants to convert an old
oil platform to an LNG terminal offshore from Oxnard; and Woodside
Energy Ltd. wants to build one near Malibu; Sempra Energy Co. has
already begun building a massive LNG terminal near Ensenada in Baja
California, Mexico.

Showboat promoter [seacoastonline.com]

While several of the details still have to be finalized, Parker has
guaranteed to put on the first Portsmouth International Boat Show from
May 31 to June 1. He expects the 400 exhibitors he has contacted to
draw upward of 40,000 people to the city.

I just got initial database support working for the noaadata/ais submodule (version 0.13). This has only passed one test case and I am sure there are a few gotchas waiting. Here is an example of its use on the command line. I received this AIS message today:

!AIVDM,1,1,,B,14QsCc0001JA7Cv>gfkhJa?@0D0w,0*71

Now I can prepare for the database:

./ais_msg_1.py --sql-create

This emits the database table creation command:

CREATE TABLE position ( \
	    MessageID INTEGER, RepeatIndicator INTEGER, \
	    UserID INTEGER, NavigationStatus INTEGER, \
	    ROT INTEGER, SOG DECIMAL(4,1), PositionAccuracy INTEGER, \
	    Position_longitude DECIMAL(8,5), \
	    Position_latitude DECIMAL(8,5), COG DECIMAL(4,1), \
	    TrueHeading INTEGER, TimeStamp INTEGER, \
	    RegionalReserved INTEGER, Spare INTEGER, RAIM BOOL, \
	    syncstate INTEGER, slotoffset INTEGER );

That is pretty basic and does not include a primary key, but that is enough to test a database. Now we need the insert string, which can be had by this command:

./ais_msg_1.py -d -T sql -t nmea \
       '!AIVDM,1,1,,B,14QsCc0001JA7Cv>gfkhJa?@0D0w,0*71'

This results in an insert string like this:

INSERT INTO position (RegionalReserved,NavigationStatus,\
  COG,SOG,slotoffset,Position_longitude,RepeatIndicator,\
  UserID,RAIM,syncstate,TimeStamp,Spare,PositionAccuracy,\
  MessageID,TrueHeading,ROT,Position_latitude) 
  VALUES \
  (0,0,10.6,0.1,10247,-80.147415,0,304010156,0,0,40,0,0,1,\
  295,0,25.77015833333333333333333333);

Now paste those two commands into an sqlite database:

sqlite3 ais.db3
[paste CREATE and then INSERT here]

Now you have a database.

sqlite> SELECT UserID FROM position;
304010156

I have been trying to figure out how to use scipy's faster arrays in C so that I can just access the data as a casted array from C (Barry's suggestion.) The documentation for this stuff is a twisted maze of headers with very few comments and no autogenerated summary docs or clear examples. Perhaps there is a web page or two that I need to discover, so if you know of a great example, please let me know! Here is what I did so far to prove the concept and try to do this correctly and compliant with the latest APIs. I could easy be wrong, but it "works for me" with python 2.5 and...

 i   numarray-py25    1.5.2-1      Multidim'l numerical arrays for Python
 i   scipy-core-py25  1.0.1-1      Replacement for Numeric python
 i   scipy-py25       1:0.5.2-3    Scientific tools for Python

First the python code to create arrays and pass them into C.

#!/usr/bin/env python
import numpy
a = numpy.core.array((1.1,2.2,3.3,4.9))
print 'len a in py:',a,type(a)
import karray
karray.try1(a)
print a
######################
import array
a2 = array.array('d',(4.5,6.7))
print 'len a2 in py:',a2,type(a2)
karray.try1(a2)
print a2

Both of the above ways work to create an array for the C code. Is one of these better than the other?

#include <assert.h>
#include <Python.h>
#include <numarray/libnumarray.h>
static PyObject*try1(PyObject* self, PyObject* args) {
    PyObject *anArray=0;
    if(!PyArg_ParseTuple(args, "O", &anArray)) return 0;
    printf ("isBuffer: %s\n",(isBuffer(anArray)?"yes":"no"));
    printf ("getBufferSize: %d\n", getBufferSize(anArray));
    {
	void *buffer=0;
	double *bufferDouble=0;
	long buffersize=-1;
	size_t len=-1;
	if ((buffersize = NA_getBufferPtrAndSize(anArray, 0, &buffer)) < 0) {
	    fprintf(stderr,"NA_getBufferPtrAndSize failed\n");
	    return 0;
	}
	len = buffersize/sizeof(double);
	printf ("buffersize: %ld -> %d\n",buffersize,(int)len);
	assert(buffer);
	bufferDouble = (double *)buffer;
	printf ("%lf\n",bufferDouble[0]);
	bufferDouble[0]=-99.3;  // Prove that I can change a value
    }
    return Py_BuildValue("i",1); // Dummy return of an int
}
static PyMethodDef karrayMethods[] = {  
    {"try1", try1, METH_VARARGS, "FIX: put doc here"},
    {0}
};
void
initkarray(void)   
{
    Py_InitModule("karray", karrayMethods); 
    import_libnumarray();
}

The above try1 function takes an array of type array.array or numpy.ndarray and then starts working with the object without much error checking. Inside the inner {}, you can see that I get the buffer and its length. Then I assume that a 1D double array is passed into the function. I set the first value in the array to show that changes can be passed back from C to Python in the array. If you use the array.array style, make sure to use 'd' for doubles otherwise you will see strange results. 'f' will not work right.

len a in py: [ 1.1  2.2  3.3  4.9] 
isBuffer: yes
getBufferSize: 32
buffersize: 32 -> 4
1.100000
[ 99.3   2.2   3.3   4.9]
========
len a2 in py: array('d', [4.5, 6.7000000000000002]) 
isBuffer: yes
getBufferSize: 16
buffersize: 16 -> 2
4.500000
array('d', [99.299999999999997, 6.7000000000000002])

This Google Earth demo by Schwehr and Weber illustrates ideas about visualizing acoustic tracking of whale calls and integration with visual overflight observations. I must emphasize that this is NOT real data. This is my first try at including a Screen Overlay in GE and a view point. The overlay was pretty easy, but the view point control took a few tries to get right.

WhaleLocalizationDemo-Jan2007.kmz

UNH had another power outage today. The story I heard is that a transformer went at the TV station and took out most of the campus with it.

Want to know what your mac's optical drive is capable of? I just found this command yesterday and here are the results from my 4xG5 desktop.

drutil info
 Vendor   Product           Rev 
 HL-DT-ST DVD-RW GWA-4165B  C006
   Interconnect: ATAPI
   SupportLevel: Apple Shipping
          Cache: 2048k
       CD-Write: -R, -RW, BUFE, CDText, Test, IndexPts, ISRC
      DVD-Write: -R, -RW, +R, +R DL, +RW, BUFE, Test
     Strategies: CD-TAO, CD-SAO, CD-Raw, DVD-DAO

PaddingTravelers on NOAA chart online.

It this really a NOAA server? EyeSpy [205.156.4.52]

Maybe it is some sort of hosting thing. But you can get there from NOAA's Atlantic Coast Chart On-Line Viewer.

traceroute 205.156.4.52
...
10  noaa-rtr.maxgigapop.net (206.196.177.98)  12.234 ms *  12.204 ms
11  140.90.75.131 (140.90.75.131)  12.345 ms  12.842 ms  12.360 ms
12  140.90.111.36 (140.90.111.36)  12.423 ms  12.355 ms *
*
...

OGLE Earth says Google has updated the KML documentation.

KML Samples

MIME types for Apache. This I need to try out.

    AddType application/vnd.google-earth.kml+xml .kml
    AddType application/vnd.google-earth.kmz .kmz

They still have this strange rsync software delivery process which makes it impossible to integrate with fink, but I do appreciate the Mac OSX support.

===================================================
ISIS 3.1.5 RELEASE 2007/01/17 - Mac OS X Compatible
===================================================
A new version of Isis 3 is now available. For download and installation 
information please see the Isis 3 Installation Guide:
.
    http://isis.astrogeology.usgs.gov/documents/InstallGuide 
.
To view the release notes, click on the following link:
.
    http://isis.astrogeology.usgs.gov/IsisSupport/viewtopic.php?t=1130
.
If you are only interested in HiRISE processing, and/or are restricted 
as to the amount of space you have available for ISIS and its ancillary 
data, then you can choose to only install the base and MRO specific 
data areas. To do this, use the following commands: 
.
    - cd $ISISROOT/../data 
    - rsync -avz isisdist.wr.usgs.gov:isis3data/data/base . 
    - rsync -avz isisdist.wr.usgs.gov:isis3data/data/mro . 
.
If you have a multi-OS environment, it is possible, and recommended, 
that you set it up so all versions share the same data area. To do this, 
simply set the environment variable "ISIS3DATA" to the equivalent of 
"$ISISROOT/../data" *after* sourcing the ISIS3 startup file.

I think AccuWeather needs to link about their algorithm for displaying the low. When the low of the day is above the current temperature, perhaps something needs to be changed

Sierra Club to show 'Inconvenient Truth' [seacoastonline.com]

Join Seacoast Sierra Club to watch "An Inconvenient Truth" -- Al Gore's
shocking and inspiring movie about global warming and what we can do
about it -- Tuesday, Jan. 16, at 7 p.m. in the Dover CIty Hall
auditorium. Everyone is welcome.

From ProgrammableWeb Mashup Detail: NOTAM Ace:

Notices to Airmen and Temporary Flight Restrictions using Google Map API [notamace.com]

Panbo has an article on the Iris PC Radar. The big questions for me are how much and does this think have linux drivers?

Construct sounds pretty exciting! Can I use this to declare the contents of an AIS message at the bit level? Will it be much faster than BitVector? From the web page...

Being declarative means you directly define the data structure, not the code that parses it. This means you can export your definition to ML, or write a code generator that takes a construct and converts it to a C module, just to name a few ideas.

Features

• Bit-level - you can work with individual bits and non-aligned fields without hassle
• Byte ordering - supports both little and big endian byte ordering
• Componentized - combine the primitives into more complex constructs
• Symmetrical - constructs can both parse (create an object from raw data) and build (create raw data from an object)
• Ready - comes with many built-in primitives and inventory constructs that you can use out-of-the-box.
• Declarative - you write a data structure, not code.
• Easy to debug - much easier to correct a data structure instead of code
• Easy to test - componentization means that once a component works, it always works. You can easily share constructs between modules without something breaking.
• Easy to extend - when the need arises, just subclass your very own construct
• Thread-safe - Construct is fully thread safe; you can parse the same construct in parallel threads.

Compare different datasets... www.flashearth.com

I am just about to check in gdal 1.4.0 into fink. Why is this so interesting? This lets anyone create a GoogleEarth S57 ENC in just a couple minutes.

For Mac users...

fink selfupdate
fink install gdal

Go to NOAA ENC. Scroll way down to PROCEED. Select the "Textual" interface. Change the "Show" drop down from 20 to 617. Search for Cape Cod. Check the box next to "Bay of Fundy to Cape Cod - US3EC10M". Press "Review Selections" at the top of the page. Press "Order Selections." Press "Download."

Now unzip the file. "cd ENC_ROOT/US3EC10M/" and run this command (once you have gdal 1.4.0 installed).

ogr2ogr -f KML US3EC10M.kml US3EC10M.000

Now run "open US3EC10M.kml" and turn off all layers accept DSID, CBLSUB, COALNE, RECARE, RESARE, TSSBND, TSSLPT, and USSARE.

http://estuarinebathymetry.noaa.gov/finddata.html - Data in Google Earth format. Stops just at the point where I need the data. But still very cool.

What is the difference between the SpaceNavigator SE and PE? The SE gets you support phone calls for an extra $50. Same hardware.

I was all excited to install a SpaceNavigator on my Mac. I had to install the deta driver (3dxwaremacPE-v0_1_0.dmg) and restart my mac. Now sketch crashes and I find this lovely set of messages in the system log.

3DxSketchup Plug-in: [74] Initializing...
3DxSketchup Plug-in: [95] 3Dconnexion device init OK? 1
tdxInit OK.  clientid=8198
2007-01-12 09:41:00.231 3DxSketchUpHelper[620] #3DxSketchUpHelper#[34]
imgPath=/Library/Application Support/Google SketchUp
6/SketchUp/Plugins/3DxSketchUp.bundle.plugin/Contents/MacOS/3DxSketchUpHelper.app/Contents/Resources/banner.jpg
3DxSketchup Plug-in: active model not available
3DxSketchup Plug-in: active model not available
3DxSketchup Plug-in: [117] Initialization complete.
free_Command
free_Command
free_Command
free_Command
#3DxSketchUpHelper#
 /Users/ep/svn_seefeld/3dxsketchup/trunk/extension/macdev/../src/mac/HelperController.m[56]
 3DxSketchUpHelper is exiting...
3DxSketchup Plug-in: active model not available
terminate called after throwing an instance of 'std::logic_error'
  what():  3DxSketchup Plug-in: active model not available

If I unplug the SpaceNavigator, start SketchUp 6, add come model elements, then plug in the SpaceNavigator, then it works. The mapping is exactly opposite of what I expect. You move the model, not the camera, so I can barely move around with the device. It takes about 10 minutes to get used to this mapping.

After that first time of removing the SpaceNavigator, I can now kill and start SkeptUp without any troubles while leaving the device plugged in.

This is all using a 4xG5 2.5 GHz PowerPC machine running 10.4.8 and Google SketchUp 6.0.277.

Check out the Husky airboat that is good for on the ice. Ice no obstacle for new rescue boat [NBC Saint Paul] How does the hull look after a few weeks of using it?

NOAA Voluntary Observing Ship (VOS) Weather Feed for Google Earth [destinsharks.com]

Wow... this is kind of how I have been feeling about SOAP. My 1st day with SOAP was great. Then I learn that the SOAP library I was using is no longer supported (SOAPpy). Then I learned that ZSI and Apache Axis do not play well together. Then I tried just creating a SOAP request and parsing the reply myself. Then I tried OPeNDAP.

The S Stands For Simple

Here is about the simplest possible C module (taken from Extending Python with C: Part 1[pyzine]. First the C source (tiny.c):

#include <Python.h>
/* the function that is callable from python */
static PyObject*
tiny(PyObject* self, PyObject* args)   /* proto for all exposed funcs */
{
    if(!PyArg_ParseTuple(args, ""))    /* receive arguments (no args) */
	return 0;                      /* propagate error if any */
    return Py_BuildValue("s","it works!"); /* build and return result */
}
/* Lookup table of functions */
static PyMethodDef tinyMethods[] = {   /* methods exposed from module */
    {"tiny", tiny, METH_VARARGS, "A tiny function"},    /* descriptor */
    {0}                                                 /* sentinel   */
};
/*  Hook into python */
void
inittiny(void)                             /* called by Python on import */
{
    Py_InitModule("tiny", tinyMethods);  /* init module with methods */
}

Here is the distutils based setup.py file that compiles the code.

#!/usr/bin/env python
from distutils.core import setup, Extension
setup(
    name="samples", 
    version="0.1",
    ext_modules = [
	Extension("tiny",["tiny.c"])
	]
    )

Now install and run the module. This is for python2.5 on fink/mac osx. The path will be different on other systems.

chmod +x setup.py
mkdir inst
./setup.py install --root=inst
  running install
  running build
  running build_ext
  building 'tiny' extension
  gcc -fno-strict-aliasing -Wno-long-double -no-cpp-precomp -mno-fused-madd -DNDEBUG -g -O3 -Wall -Wstrict-prototypes -I/sw/include/python2.5 -c tiny.c -o build/temp.macosx-10.4-ppc-2.5/tiny.o
g  cc -L/sw/lib -bundle -undefined dynamic_lookup build/temp.macosx-10.4-ppc-2.5/tiny.o -o build/lib.macosx-10.4-ppc-2.5/tiny.so
  running install_lib
  copying build/lib.macosx-10.4-ppc-2.5/tiny.so -> inst/sw/lib/python2.5/site-packages
  running install_egg_info
  Removing inst/sw/lib/python2.5/site-packages/samples-0.1-py2.5.egg-info
  Writing inst/sw/lib/python2.5/site-packages/samples-0.1-py2.5.egg-info



PYTHONPATH=inst/sw/lib/python2.5/site-packages ipython
In [1]: import tiny
In [2]: tiny.tiny()
Out[2]: 'it works!'

Tiny has limited documentation.

help (tiny.tiny)
Help on built-in function tiny in module tiny:



tiny(...)
    A tiny function

Google Earth 4 is now officially released. My copy on the mac says it is version 4.0.2722.

http://earth.google.com/download-earth.html

Available in June and definitely out of my price range. This text from their website is a little strange. Would have thought that they would have gotten approval first. Perhaps the process would have blown the secrecy.

This device has not been authorized as required by the rules of
the Federal Communications Commission. This device is not, and may not
be, offered for sale or lease, or sold or leased, until authorization
is obtained.

Integrating existing C/C++ code or converting python to C/C++ are common tasks for python folks. There are many ways to attack this problem. Which one you use depends on many factors, but here are some starting points.

• ctypes
  • ctypes at sf - pre-inclusion to python 2.5
  • ctypes in python 2.5
  • header to ctypes
• SWIG - build interfaces for many languages
• Extending and Embedding the Python Interpreter [docs.python.org]
• Python/C API Reference Manual [docs.python.org]
• SIP - created for PyQT
• PyCXX - inline C code
• Weave - another inlining package [scipy.org]
• Instant inlining of C/C++ code in Python. ("built on top of SWIG" but claims to only work on unix systems)
• PyRex
• PyInline

That is quite a list of options.

From the Fairplay Daily News Headlines:

US nuclear sub hits Japanese VLCC
NO injuries have been reported after a US Navy submarine collided with
a Japanese VLCC south of the Straits of Hormuz last night.

I haven't even had a chance to touch the last version that is installed on my machine, but Google SketchUp 6 was just released.

SpaceNavigator Works on Mac OS X [GoogleEarthBlog]

Have to grab the one that our lab has and give it a try today.

Check out these two courses at MIT... How different would these courses be if they were taught in 2007 instead of 2003/4?

Spatial Database Management and Advanced Geographic Information Systems

Environmental Engineering Applications of Geographic Information Systems

Google Earth now has texture mapped buildings. Go google! Here is on in San Francisco. Just wish they would get to the Hampton Roads/Norfolk area!

Looks like Google is still working out deployment issues with Google Earth. Yesterday, I checked for updates to GE on the mac and GE claimed there were none. Knowing otherwise, I downloaded a new version from the web. Today, I tried the same on another computer. I got the download, installed it, and it immediately said that I should update... wow. But my copy was just 2 minutes old. Then when I down loaded what was supposed to be 4.0.2723, it claims to be 4.0.2722 with a build date of a few days ago.... and it does not prompt me for a new version. When I go into check of updates, it wants me to download yet another dmg. This is quirky, but not really a big deal.

JJ has been looking into handling byte order issues between Mac OSX on PPC and x86. I only have a PPC mac right now. Here is my quick exploration of endian. I should see what works on an x86 box, but I'm not on the net, nor do I have the qemu simulator installed that would give me a whole bunch of os images to try. Is there a vms emulator out there? That would be a twisted flashback to the late 80's (e.g. {hal,gal}.arc.nasa.gov to anyone who remembers those boxes).

I dug through the headers on my system and found these macros. I don't know if any are a part of a standard.

// g++ -o endian -Wall -g3 endian.C
#include <iostream>
#include <sys/sysctl.h>
using namespace std;
int main(int argc, char *argv[]) {
    // sys/wait.h
    cout << "__DARWIN_BYTE_ORDER    " << __DARWIN_BYTE_ORDER << "\n";
    cout << "__DARWIN_LITTLE_ENDIAN " << __DARWIN_LITTLE_ENDIAN << "\n";
    cout << "__DARWIN_BIG_ENDIAN    " << __DARWIN_BIG_ENDIAN << "\n";
    // sys/sysctl.h
    cout << "HW_BYTEORDER           " << HW_BYTEORDER << "\n";
    // libkern/OSByteOrder.h
#ifdef __LITTLE_ENDIAN__    
    cout << "__LITTLE_ENDIAN__      " <<__LITTLE_ENDIAN__  << "\n";
#endif
#ifdef __BIG_ENDIAN__    
    cout << "__BIG_ENDIAN__         " << __BIG_ENDIAN__ << "\n";
#endif
    return EXIT_SUCCESS;
}

Here is what I get on my G4:

./endian 
__DARWIN_BYTE_ORDER    4321
__DARWIN_LITTLE_ENDIAN 1234
__DARWIN_BIG_ENDIAN    4321
HW_BYTEORDER           4
__BIG_ENDIAN__         1

I really wish that there was a standard unix command line tool that gave details on byte order and formats for both integers and floats.

For reference, here is a program from my density package tells the endian for a machine. It doesn't handle some of the stranger endians.

#include <cstdio>
#include <cstdlib>
using namespace std;
int main(void) {
  int s;
  char *c=(char *)&s;
  c[0]=0x00; c[1]=0x01; c[2]=0x02; c[3]=0x03;
  if (0x00010203==s) {printf("BIGENDIAN"   );return(EXIT_SUCCESS);}
  if (0x03020100==s) {printf("LITTLEENDIAN");return(EXIT_SUCCESS);}
  printf("ERROR!!!  I can not cope with this beast... HELP!!!!\n");
  return (EXIT_FAILURE);
}

There are also are all the byte order handling like utonl, htons, ntohl, and ntohs.

If someone has a really great reference on dealing with these nagging issues, please send me a link to it!

I got inspired to play with my point and shoot Nikon Coolpix 4600 by jwz's macro photos. My camera and lense are not even close to the quality of his, but it was still fun to try it out.

Schneier on Security: Licensing Boaters

There are lots of good reasons to license boats and boaters, just
as there are to license cars and drivers. But counterterrorism is not
one of them.

Saw this linked from JWZ's livejournal site. Nice to see some of the content that was on reality.sgi.com come back from the dead.

A Multifocus Method for Controlling Depth of Field page by Paul Haeberli. Similiar to the depth from focus done by Mark Lemmon with the MER microscopic imager (accept that the camera was moved, not the focus of the camera).

This Dashboard image does not even say it all. It is over 60 degrees F and sunny. Is this really New Hampshire???

I just gave the waterlevel messages a go to verify that I could take a NMEA string for a binary message and extract a kml placemark.

First get the binary payload from the msg 8 packet:

./ais_msg_8.py -d -t nmea '!AIVDM,1,1,,A,8E2HH=dtdd<=B4qI:`Zd0000?3@qF000BeH0,0*12'
BinaryData:
        MessageID:        8
        RepeatIndicator:  1
        UserID:           338040883
        Spare:            0
        BinaryData:       000000010110111000010000000\
00001000100101011010110000000001110001101101100111100\
101100101100001100001101010010000100111001011001001010\
101000100010101100000000000000000000000000001111000011\
010000111001010110000000000000000000010010101101011000\
000000

Then, I have to paste that binary string into the waterlevel decoder:

./waterlevel.py -d -t binary -T kml-full '0000000\
101101110000100000000000100010010101101011000000000111\
000110110110011110010110010110000110000110101001000010\
011100101100100101010100010001010110000000000000000000\
000000000111100001101000011100101011000000000000000000\
0010010101101011000000000' >  wl-8632200.kml
open  wl-8632200.kml

With that, I have a nice view of the tide guage data via AIS (ais-over-internet).

wl-8632200.kml

I just got kml generation working for my noaadata project. The code generator looks for the first fields names containing "longitude" and "latitude" and uses these for positioning the message pop-up. Expect to see this work in the noaadata-0.9 release. It also looks for a titlefield attribute in the message field. It then uses this subfield for the name field within the kml.

./ais_msg_1.py --decode -T kml-full -t nmea '!AIVDM,1,1,,A,15RTgt0PAso;90TKcjM8h6g208CQ,0*4A,r003679900,1168010561' > 371798000.kml

View the resulting GoogleEarth file: 371798000.kml

My new noaadata package includes ais message handling. Here is a quick demo of what it can do as of version 0.8. You can download it here: http://vislab-ccom.unh.edu/~schwehr/ais/waterlevel/downloads/ (Note: this location may change in the future).

First, here is an N-AIS ship position message in its native NMEA-AIS format. I grabbed this one at random from the stream. I don't know what ship this is or where it is located.

!AIVDM,1,1,,A,15RTgt0PAso;90TKcjM8h6g208CQ,0*4A,r003679900,1168010561

Now, I can ask ais_msg_1.py for help.

./ais_msg_1.py -h
Usage: ais_msg_1.py [options]
Options:
  --version             show program's version number and exit
  -h, --help            show this help message and exit
  --doc-test            run the documentation tests
  --unit-test           run the unit tests
  -v, --verbose         Make the test output verbose
  -t IOTYPE, --type=IOTYPE
                        What kind of string to expect (binary, nmeapayload,
                        nmea) [default: nmeapayload]
  -T OUTPUTTYPE, --output-type=OUTPUTTYPE
                        What kind of string to output (std, html, xml)
                        [default: std]
  -o OUTPUTFILENAME, --output=OUTPUTFILENAME
                        Name of the python file to write [default: stdout]
  -d, --decode          decode a "position" AIS message
  -e, --encode          encode a "position" AIS message
  --RepeatIndicator-field=uint
                        Field parameter value [default: 0]
  --UserID-field=uint   Field parameter value [default: none]
...

I want to decode the message to see what is in there:

./ais_msg_1.py -d -t nmea \
'!AIVDM,1,1,,A,15RTgt0PAso;90TKcjM8h6g208CQ,0*4A,r003679900,1168010561'



position:
        MessageID:           1
        RepeatIndicator:     0
        UserID:              371798000
        NavigationStatus:    0
        ROT:                 -127
        SOG:                 12.3
        PositionAccuracy:    1
        Position_longitude:  -123.395383333
        Position_latitude:   48.38163333333
        COG:                 224
        TrueHeading:         215
        TimeStamp:           33
        RegionalReserved:    0
        Spare:               0
        RAIM:                False
        syncstate:           0
        slotoffset:          4252

Note that slotoffset is not right. The library can't currently decode conditionals. Also, the rate of turn is not right either. The library does not support decodings other than linear. ROT requires a sqrt.

What if I want to see a pretty html for that has the entries decoded by lookup tables if possible? That is now trival to do!

./ais_msg_1.py -d -t nmea -T html '!AIVDM,1,1,,A,15RTgt0PAso;90TKcjM8h6g208CQ,0*4A,r003679900,1168010561' > msg.html
open msg.html

Here is what msg.html contains. I did a a couple quick edits to make it look nicer in the blog. First, I removed all extra white space in the html. Second, I shorted the decimal values. Third, I changed a header 3 tag to a bold tag.

position

Field Name	Type	Value	Value in Lookup Table
MessageID	uint	1	1
RepeatIndicator	uint	0	default
UserID	uint	371798000	371798000
NavigationStatus	uint	0	under way using engine
ROT	int	-127	-127
SOG	udecimal	12.3	Missing LUT entry
PositionAccuracy	uint	1	high (greater than 10 m)
Position_longitude	decimal	-123.3953833	-123.3953833
Position_latitude	decimal	48.38163333	48.38163333
COG	udecimal	224	224
TrueHeading	uint	215	215
TimeStamp	uint	33	Missing LUT entry
RegionalReserved	uint	0	0
Spare	uint	0	0
RAIM	bool	False	not in use
syncstate	uint	0	UTC direct
slotoffset	uint	4252	4252

And finally, we should have Google Earth show us where that was.

Finding this made my day. Noaa has NOAA RSS which should make an easy way to stay connected with NOAA. There are also podcasts, which I have yet to check out.

Thanks to Roberto De Almeida for replying to my posts on the pydap google group (http://groups.google.com/group/pydap), I now have a working test script that pulls just one row of data from the NOAA CO-OPS server. I have a serious speed problem (22 seconds to fetch), but this is great and I am now not depending on the obsolete SOAPpy package. Note that someone has picked SOAPpy back up and started maintaining it,but most of the developers are on to ZSI.

Here is my small test script to fetch one waterlevel point.

#!/usr/bin/env python
import dap.client
import urllib
dataset=dap.client.open('http://opendap.co-ops.nos.noaa.gov/dods/IOOS/Raw_Water_Level')
seq = dataset['WATERLEVEL_RAW_PX']
reqStr=urllib.quote('_STATION_ID="1615680"&_BEGIN_DATE="20060101 10:06"&_END_DATE="20060101 10:06"&_DATUM="MLLW"')
print 'reqStr:',reqStr
filt_seq=seq.filter(reqStr) 
data = filt_seq._get_data() # Is this bad - using private method?
print 'Found this many waterlevel points:',len(data)
print data[-1]

Many thanks to Rob Cermak on the pydap mailing list for helping me figure this out. Clearly, I have more to learn about how to do these queries efficiently. However, this is good enough for now.

#!/usr/bin/env python
'''
Try out pydap to fetch the current water level from the NOAA Co-ops
server.  This is a really slow way to pull one data point!
'''
import dap.client
dataset=dap.client.open('http://opendap.co-ops.nos.noaa.gov/dods/IOOS/Raw_Water_Level')
#
print 'Dataset keys:'
for item in dataset.keys(): print '  ',item
seq = dataset['WATERLEVEL_RAW_PX']
filt_seq=seq.filter('_STATION_ID="1615680"&_BEGIN_DATE="20060101"&_END_DATE="20060101"&_DATUM="MLLW"') 
print 'filter keys:'
for item in filt_seq.keys(): print '  ',item
#
# Print the results.  To just get field, do this...
#
print 'Just the WL_VALUE field:',filt_seq['WL_VALUE'][:][-1]
#
# Fetch all the fields...
print 'One data point:'
for item in filt_seq.keys():
    print '  ',item,':',filt_seq[item][:][-1] 

When I run the above, I get this after a long wait. I really am fetching too many points and I am refetching for each field.

./test_pydap.py 
Dataset keys:
   WATERLEVEL_RAW_PX
filter keys:
   _STATION_ID
   _DATUM
   _BEGIN_DATE
   _END_DATE
   DCP
   SENSOR_ID
   DATE_TIME
   WL_VALUE
   SIGMA
   O
   F
   R
   L
Just the WL_VALUE field: 0.474
One data point:
   _STATION_ID : 1615680
   _DATUM : MLLW
   _BEGIN_DATE : 20060101
   _END_DATE : 20060101
   DCP : 1
   SENSOR_ID : A1
   DATE_TIME : l
   WL_VALUE : 0.474
   SIGMA : 0.064
   O : 0
   F : 0
   R : 0
   L : 0

The answer on speed is to both constrain the time window to the 6 minute interval (thereby not pulling the whole day) and getting all the fields in one call. I was able to do the later from ipython like this:

filt_seq._get_data()
...
 ('1615680',
  'MLLW',
  '20060101',
  '20060101',
  '1',
  'A1',
  'Jan  1 2006 11:48PM',
  0.43099999999999999,
  0.097000000000000003,
  0,
  0,
  0,
  0),
 ('1615680',
  'MLLW',
  '20060101',
  '20060101',
  '1',
  'A1',
  'Jan  1 2006 11:54PM',
  0.47399999999999998,
  0.064000000000000001,
  0,
  0,
  0,
  0)]

For some reason, I can't do this from a script? This happens if I tack a _get_data() on the end of the above script. I get this loveliness:

Traceback (most recent call last):
  File "./test_pydap.py", line 25, in 
    data = filt_seq._get_data()
  File "/sw/lib/python2.5/site-packages/dap/dtypes.py", line 223, in _get_data
    return _build_data(self.level, *[var.data for var in self.values()])
  File "/sw/lib/python2.5/site-packages/dap/dtypes.py", line 517, in _build_data
    out = [_build_data(level-1, *els) for els in itertools.izip(*vars_)]
  File "/sw/lib/python2.5/site-packages/dap/proxy.py", line 37, in __iter__
    return iter(self[:])
  File "/sw/lib/python2.5/site-packages/dap/proxy.py", line 105, in __getitem__
    resp, data = openurl(url, self.cache, self.username, self.password)
  File "/sw/lib/python2.5/site-packages/dap/util/http.py", line 21, in openurl
    raise ClientError(msg)
dap.exceptions.ClientError: 'Server error 0: "java.sql.SQLException:
JZ006: Caught IOException: java.io.IOException: JZ0PA: The query has
been cancelled and the response discarded.  The cancel was probably
issued by another statement on the connection.";\n};\nError {\n code =
0;\n message = "java.sql.SQLException: JZ0S2: Statement object has
already been closed."'

I am not sure what that java SQL error is about. But when I use a very limited script, it all works?

#!/usr/bin/env python
import dap.client
dataset=dap.client.open('http://opendap.co-ops.nos.noaa.gov/dods/IOOS/Raw_Water_Level')
seq = dataset['WATERLEVEL_RAW_PX']
filt_seq=seq.filter('_STATION_ID="1615680"&_BEGIN_DATE="20060101"&_END_DATE="20060101"&_DATUM="MLLW"') 
data = filt_seq._get_data()
print 'Found this many waterlevel points:',len(data)
print data[-1]

Results in this:

./test_pydap2.py 
Dataset keys:
   WATERLEVEL_RAW_PX
Found this many waterlevel points: 240
('1615680', 'MLLW', '20060101', '20060101', '1', 'A1', 'Jan  1 2006 11:54PM', 0.47399999999999998, 0.064000000000000001, 0, 0, 0, 0)

Now to jump back to that performance problem. I only want one sample, not 240. Pulling only one should be much faster. If I try limiting the start and end dates by adding time, that should work, but it doesn't.

filt_seq=seq.filter('_STATION_ID="1615680"&_BEGIN_DATE="20060101 10:00"&_END_DATE="20060101 11:00"&_DATUM="MLLW"') 

Gives another server side error.

dap.exceptions.ClientError: 'Server error 3: "Encountered
\'\\\'20060101\' at line 1, column 65.\nWas expecting one of:\n \'*\'
...\n \'{\' ...\n  ...\n  ...\n  ...\n  ...\n
";\n};\nError {\n code = 0;\n message = "java.sql.SQLException: JZ0S2:
Statement object has already been closed."'

When formatted for humans...

dap.exceptions.ClientError: 
Server error 3: 
"Encountered '\'20060101\' at line 1, column 65.
Was expecting one of:
    \'*\' ...
    \'{\' ...
     ...
     ...
     ...
     ...
    ";
};
Error {
    code = 0;
    message = "java.sql.SQLException: JZ0S2: Statement object has already been closed."

Does this mean there is an error in the server side parser? The spec (Raw_Water_Level.das) says that this should be okay, right?

_BEGIN_DATE {
        String long_name "REQUIRED 8- to 14-character Begin Date (yyyymmdd hh:mi) surrounded in double quotes";
        String Examples " 20040321  OR 20030712 19:02 ";
    }

I do have a 14 character string. Hmmm. Rob Cermak had an example URL for pulling the data. Time to try out the URL query method. Try it in a pop-up window.

http://opendap.co-ops.nos.noaa.gov/dods/IOOS/Raw_Water_Level.ascii?\
WATERLEVEL_RAW_PX.WL_VALUE&WATERLEVEL_RAW_PX._STATION_ID="1615680"&\
WATERLEVEL_RAW_PX._DATUM="MLLW"&WATERLEVEL_RAW_PX._BEGIN_DATE="20060101"&\
WATERLEVEL_RAW_PX._END_DATE="20060102"&WATERLEVEL_RAW_PX.WL_VALUE!=-999.9

You can try that on the command line something like this, but remove the \ characters and make it one monster line.

wget -O data 'http://opendap.co-ops.nos.noaa.gov/dods/IOOS/Raw_Water_Level.ascii?\
WATERLEVEL_RAW_PX.WL_VALUE&WATERLEVEL_RAW_PX._STATION_ID="1615680"&\
WATERLEVEL_RAW_PX._DATUM="MLLW"&\
WATERLEVEL_RAW_PX._BEGIN_DATE="20060101"&\
WATERLEVEL_RAW_PX._END_DATE="20060102"&\
WATERLEVEL_RAW_PX.WL_VALUE!=-999.9'
# Then check out the results...
wc -l data 
488 data
head data
Dataset {
    Sequence {
        Float64 WL_VALUE;
    } WATERLEVEL_RAW_PX;
} WATERLEVEL_RAW_PX;
=============================================
WATERLEVEL_RAW_PX.WL_VALUE
0.436
0.466
0.428

Note that I switched out the - for = to not freak out nano blogger. Now I should be able to specify time and get fewer results.

wget -O data.1hr 'http://opendap.co-ops.nos.noaa.gov/dods/IOOS/Raw_Water_Level.ascii?\
WATERLEVEL_RAW_PX.WL_VALUE&
WATERLEVEL_RAW_PX._STATION_ID="1615680"&\
WATERLEVEL_RAW_PX._DATUM="MLLW"&\
WATERLEVEL_RAW_PX._BEGIN_DATE="20060101"&\
WATERLEVEL_RAW_PX._END_DATE="20060102"&\
WATERLEVEL_RAW_PX.WL_VALUE!=-999.9'
# Then check out the results...
wc -l data.1hr
19 data.1hr

I can also pull the times for those samples.

wget -O data.1hr 'http://opendap.co-ops.nos.noaa.gov/dods/IOOS/Raw_Water_Level.ascii?\
WATERLEVEL_RAW_PX.DATE_TIME&
WATERLEVEL_RAW_PX._STATION_ID="1615680"&\
WATERLEVEL_RAW_PX._DATUM="MLLW"&\
WATERLEVEL_RAW_PX._BEGIN_DATE="20060101"&\
WATERLEVEL_RAW_PX._END_DATE="20060102"&\
WATERLEVEL_RAW_PX.WL_VALUE!=-999.9'

Which give:

Dataset {
    Sequence {
        String DATE_TIME;
    } WATERLEVEL_RAW_PX;
} WATERLEVEL_RAW_PX;
===================================
WATERLEVEL_RAW_PX.DATE_TIME
"Jan  1 2006 10:00AM"
"Jan  1 2006 10:06AM"
"Jan  1 2006 10:12AM"
"Jan  1 2006 10:18AM"
"Jan  1 2006 10:24AM"
"Jan  1 2006 10:30AM"
"Jan  1 2006 10:36AM"
"Jan  1 2006 10:42AM"
"Jan  1 2006 10:48AM"
"Jan  1 2006 10:54AM"
"Jan  1 2006 11:00AM"

I do not understand why I can't get a smaller time window through dap, but it does work through a URL.

Have I lost you yet???