On 26-Apr-2012, Sibyll Schaphoff wrote: also die Skripts inklusive Daten findet ihr unter: /iplex/01/elis/sibyll/hwsd wobei process_hwsd.r die Original-HWSD vorverarbeitet (von 30 sec auf 0.5 Grad) und create_soil_map.r den LPJ input erzeugt. Fuer die Funktion TT.points.in.classes in diesem Skript muss das R Paket "soiltexture" geladen werden. Doku gibts hier: http://cran.r-project.org/web/packages/soiltexture/soiltexture.pdf Stefan Petri: Dank Visualisierungsserver mit genug Hauptspeicher hab ich die Daten jetzt auch in NetCDF gewandelt bekommen. # #unzip # #BYTEORDER I #LAYOUT BIL #NROWS 21600 #NCOLS 43200 #NBANDS 1 #NBITS 16 #BANDROWBYTES 86400 #TOTALROWBYTES 86400 #BANDGAPBYTES 0 # # # 0.00833333333333 # 0.00000000000000 # 0.00000000000000 # -0.00833333333333 # -179.99583333333334 # 89.99583333326137 # #1/3600*30=.00833333333333333310 # #set memory/size=15000 #define axis/x=179.99583333333334w:179.99583333333334e:0.00833333333333 x30sec #define axis/y=89.99583333326137S:89.99583333326137N:0.00833333333333 y30sec_1 #define grid/x=x30sec/y=y30sec g30sec #file/var=hwsd_rev/grid=g30sec/format=stream/type=i2 "hwsd.bil" ##sigh. This data is upside down so far. Use samplej() to reverse the lat axis. ##ncpqd -a -lat reverses both, axis and data, thus no real improvement. # Using the /order parameter with above file command did not work. #let reverse_j = 21600-j[g=g30sec] #let hwsd=samplej(hwsd_rev,reverse_j) #set variable/title="Harmonized World Soil Database (version 1.21)" hwsd #save/file="/scratch/01/petri/hwsd.cdf" hwsd # Andere Moeglichkeit zur Visualisierung: rawtopgm -maxval 65535 -bpp 2 43200 21600 < hwsd.bil > hwsd.pgm erzeugt ein 16-bit Portable Gray Map. Das kann angeschtu werden mit xv, wenn auch nur im 8-bit-colormap-Modus. Imagemagick (display) bricht mit error ab. GIMP als poor man's GIS missbrauchen klappt leider auch nicht. GIMP bis version 2.6 kann nur 8-bit PGMs lesen. Evtl. kann GIMP version 2.8 das, die ist aber am PIK nirgends (?) installiert. Und die zugehoerigen Metadaten im MS-Access-Format lassen sich auf dem wincenter bearbeiten und nach .csv-Format exportieren (HWSD.mdb -> HWSD_DATA.csv). Zusammen mit den R-Scripten werd ich jetzt also etwas bauen, das passende LPJ-Input-Dateien mit land_lad-kompatibler Landmaske erstellt. Das Koordinaten-File grid.bin habe ich schon. Das hat 90881 Land-Zellen, inklusive Antarktis etc. Wenn ich das richtig verstehe, will create_soil_map.r eine schon vorhandene soil.bin einelesen, und die dann modifizieren anhand der HWSD-Daten? Aber ich habe ja eigentlich keine passende alte soil.bin zur Hand, sondern wuerde eigentlich ``einfach'' eine neue from scratch bauen. Und ich sehe dort vector.icerock[i] <- frac.RK[ilon[i],ilat[i]] + frac.GG[ilon[i],ilat[i]] Aber fuer die Berechnung der Albedos (und vermutlich spaeter auch andere Sachen) haette ich gerne Eis-Zellen und Felsen-Zellen getrennt aufgefuehrt. Ich habe erst das process_hwsd.r laufen lassen, und jetzt das create_soil_map.r . #handle borderline cases soil.class[which(soil.class[,"Lo"] & soil.class[,"SiLo"]),"SiLo"] <- FALSE soil.class[which(soil.class[,"SaClLo"] & soil.class[,"SaLo"]),"SaLo"] <- FALSE soil.class[which(soil.class[,"Cl"] & soil.class[,"ClLo"]),"Cl"] <- FALSE soil.class[which(soil.class[,"SaLo"] & soil.class[,"LoSa"]),"SaLo"] <- FALSE Da ist mir nun vollkommen unklar, wie ich, statt eine alte soil.bin zu modifizieren, eine komplett neue bauen kann. Mit image.plot() kann ich mir die data.XXX und frac.XXX die aus dem HWSD-File hergeleitet sind, anschauen. Probleme dabei: - die Umrisse haben nichts mit der Landmaske im eingelesenen grid.bin zu tun - ich finde nirgends Daten fuer die Antarktis. frac.GG enthaelt nur Groenland und etwas Kleckerkram in Bergregionen. Sibyll: richtig ist, dass wir keinen Soil-code fuer die Antarktis haben, und in Groenland auch nichts, weil da halt nur Eis ist und kein Boden. Sibyll: Die alte soil map wurde nur eingelesen um bessere Werte fuer nicht-uebereinstimmende Land-Masken-Punkte zu finden. Das kann fuer uns einfach raus. Das o.a. erzeugte hwsd.cdf mit ncview beschaut zeigt, dass fuer die Antarktis gar keine Bodendaten in hwsd.bil enthalten sind (alle Werte sind 0, wie fuer Ozean-Zellen). Das Groenlaendische Inland dagegen hat durchgehend Bodentyp Nr. 6998. Das ist in den Metadaten HWSD_DATA.csv ein Typ ohne Eigenschaften (6998 = MU_GLOBAL, zweite Spalte in HWSD_DATA.csv , nicht erste Spalte ID), nur mit der Kennung ``GG'' (Glacier and Permanent Snow) versehen in Spalte SU_SYM74. Sibyll: Ja genau der Bodentyp ist 6998, aber nicht die ID sondern MU_GLOBAL ist die Referenz und der hat keine Eigenschaften, wenn Du das jetzt im process_hwsd.r abfaengst und frac.GG auf 100 bzw 1 setzt sollte das Ganze gegessen sein. Im 720x360-Gitter faellt Zelle ilon=275,ilat=42 ins Groenlandeis. hwsd.bil hat fuer alle zugehoerigen Zellen (lon=16440:16500,lat=2460:2520) den Wert 6998. Loesung: im create_soil_map.r war bei der Trennung zwischen Ice und Rock ein Dreher drin. Klein aber gemein. Naechster Schritt: Lake Fraction. frac.WR in process_hwsd.r enthaelt den Anteil der Wasserflaechen. Diese in einen Cell-Vector gewandelt und in eine Datei geschrieben, analog zur Erzeugung der soil map, ergibt schon ein brauchbar aussehendes Ergebnis. Offenbar sind darin aber Rand/Binnen-Meere, die aus Sicht des Ozean-Modells zubetoniert sind (Aralsee, Rotes Meer, Golf von Aden, Persischer Golf, Great Lakes, Aermelkanal), mit Hilfe der Nachbarschaftssuche mit Soil-Werten befuellt worden? Im hwsd.bil stehen an den Stellen 0-Werte (missing data bzw. Ozean). Loesungsansatz: fuer die Binnenmeere die 0-Werte durch neu erzeugte ``Boden''-Typen ersetzen. HWSD_DATA field SU_SYM90 has alphabetically last entry named WR. HWSD_DATA field MU_GLOBAL has highest entry 32050 HWSD_DATA last field has ID 48148 -> add artificial entries to denote hand-filled continental seas. ID,MU_GLOBAL,MU_SOURCE1,MU_SOURCE2,ISSOIL,SHARE,SEQ,SU_SYM74,SU_CODE74,SU_SYM85,SU_CODE85,SU_SYM90,SU_CODE90,...,ADD_PROP,... # Great Lakes 48149,32100,,,0,100.00,,,,,,CSGL,,,,,,,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, # Aral Sea 48150,32101,,,0,100.00,,,,,,CSAS,,,,,,,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, # Red Sea 48151,32102,,,0,100.00,,,,,,CSRS,,,,,,,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, # Gulf of Aden 48152,32103,,,0,100.00,,,,,,CSGA,,,,,,,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, # Persian Gulf 48153,32104,,,0,100.00,,,,,,CSPG,,,,,,,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, # English Channel 48154,32105,,,0,100.00,,,,,,CSEC,,,,,,,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, # Denmark Street 48155,32106,,,0,100.00,,,,,,CSDS,,,,,,,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, # Gibraltar 48156,32107,,,0,100.00,,,,,,CSGI,,,,,,,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,