DZone Snippets: gis code

Download coordinate from maproom

lindenb (Pierre) — Thu, 31 May 2007 07:29:53 GMT

// description of your code here



IFS="

"

for C in Africa Antarctica Asia Australasia Europe North_America South_America

do

	mkdir -p GIS/$C

	wget -O jeter.html -q "http://www.maproom.psu.edu/cgi-bin/dcw/dcwarea.cgi?$C"

	for P in `egrep '^' -f2`

	do

		curl -d "area=$C" -d "country=$P" -o jeter2.html "http://www.maproom.psu.edu/cgi-bin/dcw/dcwcountry0.cgi"

		P2=`grep FORM jeter2.html | grep point10 | tr "<>" "\n\n" | grep count| cut -d ' ' -f4 | cut -d '=' -f2`

		#curl -d "country=$P2" -o jeter3.txt "http://www.maproom.psu.edu/cgi-bin/dcw/point10.cgi"

		wget -O jeter2.html "http://www.maproom.psu.edu/cgi-bin/dcw/point10.cgi?country=$P2"

		P3=`cat  jeter2.html | tr " >" "\n\n" | grep ftp`

		echo "##### $P2 $P3"	

		wget -O jeter3.txt "$P3"

		mv jeter3.txt > GIS/$C/${P2}.js

		sleep 3

	done

done

Location queries for a model

ajturner (Andrew Turner) — Wed, 12 Apr 2006 05:46:45 GMT

// These are class functions which are useful for finding locations in a database
// when PostGIS is not available (when using MySQL for example)



def self.deg2rad(deg)

	(deg * Math::PI / 180)

end



def self.rad2deg(rad)

	(rad * 180 / Math::PI)

end



def self.acos(rad)

	Math.atan2(Math.sqrt(1 - rad**2), rad)

end



def self.distance_in_miles(loc1, loc2)

	lat1 = loc1.latitude

	lon1 = loc1.longitude

	lat2 = loc2.latitude

	lon2 = loc2.longitude

	theta = lon1 - lon2

	

	dist = Math.sin(self.deg2rad(lat1)) * Math.sin(deg2rad(lat2)) 

		+ Math.cos(self.deg2rad(lat1)) * Math.cos(self.deg2rad(lat2)) * Math.cos(deg2rad(theta))

	

	dist = self.rad2deg(self.acos(dist))

	

	(dist * 60 * 1.1515).round #distance in miles

end



def miles_to(location)

	Location.distance_in_miles(self, location)

end



def self.locationArea(location, miles)

	radius = miles.to_f

	latR = radius / ((6076 / 5280) * 60)

	lonR = radius / (((Math.cos(location.latitude * Math::PI / 180) * 6076) / 5280) * 60)

	

	{

		:min_latitude => location.latitude - latR,

		:min_longitude => location.longitude - lonR,

		:max_latitude => location.latitude + latR,

		:max_longitude => location.longitude + lonR

	}

end



def self.location_ids_in_range(location, miles)

	la = Location.locationArea(location, miles)

	Location.find(:all,

		:select => "locations.id",

		:conditions => ["locations.latitude <= ? and locations.latitude >= ? " + \

			" AND locations.longitude >= ? and locations.longitude <= ?",

			la[:max_latitude], la[:min_latitude], la[:min_longitude], la[:max_longitude]

		]

	).collect { |l| l.id }

end



def self.locations_in_range(location, miles)

	la = Location.locationArea(location, miles)

	Location.find(:all,

		:conditions => ["locations.latitude <= ? and locations.latitude >= ? " + \

			" AND locations.longitude >= ? and locations.longitude <= ?",

			la[:max_latitude], la[:min_latitude], la[:min_longitude], la[:max_longitude]

		]

	)

end

Shapefile2KML

jgoodall (Jon) — Wed, 16 Nov 2005 05:39:02 GMT



#requires and ArcGIS 9.x licenses

#import standard libraries

import sys, os, win32com.client



#import the kml library

import pyKML



#create gp object

gp = win32com.client.Dispatch("esriGeoprocessing.GPDispatch.1")



#get the input layer path

layer = "C:/myProjects/CUAHSI/WebServices/NCDC/asos-20051110/asos-20051110.shp" #sys.argv[1]



#get the lat, lon, name, and id fields

lat_field = "LATITUDE"

lon_field = "LONGITUDE"

elev_field = "STN_ELEV"

name_field = "LOCATION"

id_field = "CALLSIGN"



#create a cursor on this layer 

feats = gp.SearchCursor(layer)

feat = feats.Next()



#open a kml file to write points into

kml = pyKML.KML_File("C:/myProjects/CUAHSI/WebServices/NCDC/asos-20051110/test.kml")

kml.open_folder("Weather_Stations")



#loop through features

while feat:



    shape1 = feat.shape

    part = shape1.GetPart(0)

    

    #get he features attributes

    lat =  part.y

    lon = part.x

    elev = feat.GetValue(elev_field)

    name = feat.GetValue(name_field)

    id = feat.GetValue(id_field)



    #create a description from the name and id attributes                         

    description = name + "http://www.nws.noaa.gov/asos/"



    #create a kml placemaker object    

    kml.add_placemarker(lat,lon,0.0,description,id)

                           

    feat = feats.Next()



#close the folder

kml.close_folder()



#close file

kml.close()