# WAURISA Workshop 2010

# ogrinfo (vector data)

# list shape file attributes
ogrinfo –so counties2008.shp counties2008
 
# list PostGIS table attributes
ogrinfo –so PG:"host=127.0.0.1 user=postgres password=postgres dbname=workshop port=5432" counties –summary
 
#list personal geodatabase table attributes
ogrinfo -so streamnet_fishdist.mdb Fish_AllSpeciesCombined



# ogr2ogr (vector data

# convert shape file to KML
ogr2ogr -f "KML" newcounties.kml counties.shp

# select from shape file and write to new shape file
ogr2ogr -sql "SELECT *  FROM uninhabited WHERE AREA > 500000000" biguninhabited.shp uninhabited.shp



# gdalinfo, gdal_translate (raster data)
# list raster file attributes
gdalinfo wa_shade_1km.tif

#convert format
gdal_translate -of "png" wa_shade_1km.tif wa_shade_1km.png

# gdalwarp (raster data)
# reproject to geographic 4326
gdalwarp -t_srs "epsg:4326" wa_shade_1km.tif wa_shade_1km_4326.tiff

# cut to bounding box
gdal_translate -projwin 1372974.024716 545593.139169 2304119.628808 -7274.563261 wa_shade_1km.tif wa_shade_1km_cut.tif

# cut with polygon
gdalwarp  -of GTiff -r lanczos -cutline poly.shp hs1.tif hs_cut.tif

# generate Hillshade
gdaldem hillshade 


################################## 
# PostGIS

# Importing data into PostGIS
shp2pgsql -I -s 2285 counties2008.shp counties_pg > counties.sql
psql -U postgres -d weave –f counties.sql

# Can combine both with “|”:
shp2pgsql -I -s 2285 counties2008.shp counties_pg | psql -U postgres -d weave

# Other input formats with ogr2ogr
ogr2ogr -f "PostgreSQL" PG:"host=localhost user=postgres port=5432 dbname=workshop password=postgres" streamnet_fishdist.mdb -lco GEOMETRY_NAME=the_geom -t_srs "EPSG:2285" -nln "Fish_AllSpeciesCombined"

# plsql
psql -U postgres -d workshop

plsql commands
\dt
\d us_counties


# Simple spatial operations
#  Distance

select distance(setsrid((MakePoint(1622794, 150532)),2285),
setsrid((MakePoint(1622845, 150937)),2285));

#Output format
select askml(the_geom) from counties where name ilike 'king';
 
# Human readable geometry
select astext(the_geom) from counties where name ilike 'King';

# Transform
select astext(transform((setsrid(MakePoint(1622794, 150532),2285)),4326));

# Buffer
Select st_buffer(ST_Simplify(the_geom, 700), 9000) from counties where name ilike 'King';
# show in webmap


# Intersect
select name from counties where counties.the_geom && (setsrid((MakePoint(1622794, 150532)),2285)) 
and intersects (counties.the_geom,setsrid((MakePoint(1622794, 150532)),2285));

# Aggregate functions - Union of polygons
# Union all counties of the county polygon data set "us_counties " to create one polygon encompassing the area of the entire US.
select st_union (the_geom)
into us_border
from us_counties 


# This operation unions all individual datasets and groups them by states.
select st_union (the_geom), state_name
into us_states
from us_counties
group by state_name;


 ###############################################################
 # Project: MCI region
 # Calculate Landcover Area for all Counties in a 6 States region (and report by census FIPS code) then calculate by State

1. Input US Counties layer 
2. Shape file with coverage of imagery
3. NASA MODIS Landcover Imagery - MCD12Q1

convert county shape file to Albers Equal Area projection EPSG:102003
Note ->  batch files
ogr2ogr -t_srs "EPSG:102003" 102003_citynames.shp citynames.shp         

ogr2ogr -t_srs "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" 102003_modis_sinusoidal_grid_world.shp modis_sinusoidal_grid_world.shp 

# download multiple MODIS images 500 m landcover data 
# HDF Hierarchical Data Format
# SINOSIDAL projection

# gdal_translate to tif
# CHECK images
gdalinfo MCD12Q1.A2007001.h10v04.005.2009349191735.hdf

# subdatasetr band 1 - look for raster cell size...
gdalinfo HDF4_EOS:EOS_GRID:"MCD12Q1.A2007001.h09v05.005.2009349191722.hdf":MOD12Q1:Land_Cover_Type_1

# convert to tiff
gdal_translate -of Gtiff -a_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" HDF4_EOS:EOS_GRID:"MCD12Q1.A2007001.h09v05.005.2009349191722.hdf":MOD12Q1:Land_Cover_Type_1 new_2007_h09v05_sin.tif

# convert adn reproject to USA_Contiguous_Albers_Equal_Area_Conic
gdalwarp -s_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" -t_srs "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs " 2007_h09v05_sin.tif new_2007_h09v05_esri102003.tif

# check output
gdalinfo new_2007_h09v05_esri102003.tif

# gdalwarp to Albers proj - Note force cell size 463 m
# force cellsize 
 -tr 463.312716527916510 463.312716527916680
 
 gdalwarp -s_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" -t_srs "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"  -tr 463.312716527916510 463.312716527916680 2007_h09v05_sin.tif tr_new_2007_h09v05_esri102003.tif

# check output
gdalinfo tr_new_2007_h09v05_esri102003.tif 

# load into gvSIG

# One as an example
# mosaic to one tif all reprojected tifs
gdalwarp -s_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" -t_srs "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs " -tr 463.312716527916510 463.312716527916680 2008_*102003.tif new_tr_2008_mci_esri102003.tif

# directly fromm sinosidal
gdalwarp -s_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" -t_srs "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs " -tr 463.312716527916510 463.312716527916680 2008_*_sin.tif new_tr_2008_mci_esri102003.tif

# Check results in gvSIG

#reproject county shapefile
ogr2ogr -t_srs "EPSG:102003" us_counties102003.shp us_counties.shp

#Check
ogrinfo -so us_counties102003.shp us_counties102003

# several steps to run the araeas calcualtion
1. add vector and raster layer
2. rasterize counties layer ()
3. use counties and landcover raster to Tabulate Area
   using the sextante model builder in gvSIG

In Model builder:
Rasterization and Interpolation> Rasterize Vector Layer
Raster Categories Analysis > Tabulate Area 


Note  about OSGEO4W

and other shapefile manipulation utilities:
Mapserver utilities http://www.mapserver.org/utilities/
Shapelib tools http://shapelib.maptools.org/