Thanks again, but with gdalinfo and gdal_translate giving NULL after running to the end, i just cant do that..
2016-02-23 15:36 GMT-07:00 Chris Reudenbach <reudenb...@uni-marburg.de>: > Agus > > sorry I did miss the crucial part: > > if you are doing as suggestest you have to define manually the Lambertian > ymin xmax ymin ymax using the header information of the nc file. > > here an example for the unstaggered data > > library(gdalUtils) > library(raster) > library(proj4) > library(mapview) > > r<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif', > of="GTiff", > ot="Byte", > output_Raster=TRUE, > verbose=TRUE) > > finfo <- gdalinfo("geo_em.d01.nc") > ##extract parameters > lat_1=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#TRUELAT1", > finfo))],"=")[[1]][2]) > lat_2=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#TRUELAT2", > finfo))],"=")[[1]][2]) > lat0=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#CEN_LAT", > finfo))],"=")[[1]][2]) > lon0=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#CEN_LON", > finfo))],"=")[[1]][2]) > dx=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#DX", > finfo))],"=")[[1]][2]) > dy=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#DY", > finfo))],"=")[[1]][2]) > y=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#SOUTH-NORTH_PATCH_END_UNSTAG", > finfo))],"=")[[1]][2]) > x=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#WEST-EAST_PATCH_END_UNSTAG", > finfo))],"=")[[1]][2]) > x_0=0 > y_0=0 > > ## generate compliant prj4 string > projLcc=paste("+proj=lcc +lat_1=",lat_1," +lat_2=",lat_2, > " +lat_0=",lat0," +lon_0=",lon0, " +x_0=",x_0," > +y_0=",y_0,sep="") > > ## project centre coordinates > tr <- ptransform(cbind(lon0, lat0)/180*pi,'+proj=longlat +datum=WGS84 > +no_defs',proj) > > ## calculate extent using the Lambertian units (m) > xmin=as.integer(tr[1,1]-(x/2*dx-dx)) > xmax=as.integer(tr[1,1]+(x/2*dx-dx)) > ymin=as.integer(tr[1,2]-(y/2*dy-dy)) > ymax=as.integer(tr[1,2]+(y/2*dy-dy)) > wrfLccExt<-extent(xmin,xmax,ymin,ymax) > extent(r) <- extent(wrfLccExt) > projection(r) <- projLcc > mapview(r) > > > cheers Chris > > > > Am 23.02.2016 um 22:22 schrieb Chris Reudenbach: > >> Augustin >> >> just quick and dirty if you run gdalinfo("geo_em.d01.nc") your are >> getting the information about the corner coordinates the subdatasets and so >> on. Together with Dominiks suggestion you can do something like this: >> >> library(gdalUtils) >> library(raster) >> Sys.setenv(GDAL_NETCDF_BOTTOMUP="YES") >> wrffake<- "+proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m >> +no_defs" >> x<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif', >> of="GTiff", >> ot="Byte", >> output_Raster=TRUE, >> verbose=TRUE) >> wrfExt<-extent(-151.29639,-48.703613,12.355667,50.26619) >> extent(x) <- extent(wrfExt) >> projection(x) <- wrffake >> plot(x) >> >> Some remarks: >> (1) I just took the first pair of coordinates as derived from gdalinfo(" >> geo_em.d01.nc") >> you will find 4 different coordinate pairs (i did not proof which one is >> right >> >> The data is staggered (as outlined by Dominik) So some of the corner >> coordinates belongs to the staggered data and the others coordinates to >> the unstaggered ones. You will find them marked >> >> If you have installed the netcdf libs you easily can use ncview >> geo_em.d01.nc or ncdump -h geo_em.d01.nc to view the data or get more >> information of the header. >> >> Hope this helps >> >> cheers Chris >> >> >> >> Am 23.02.2016 um 21:11 schrieb Agus Camacho: >> >>> Thanks for that Dominik, >>> >>> Giving that projection to either the locations, the raster layer >>> generated >>> from the .nc file, or both, still did not work. I keep having locations >>> that should be on land falling far on the sea. Might this be a problem >>> derived from using raster with a file whose original grid distances are >>> not >>> constant? >>> >>> Here is a link with the original file which has the original coordinate >>> data. >>> >>> https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0 >>> >>> >>> 2016-02-23 12:07 GMT-07:00 Dominik Schneider < >>> dominik.schnei...@colorado.edu >>> >>>> : >>>> This looks like WRF <http://www.wrf-model.org/index.php> data. I just >>>> dealt with this. >>>> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere >>>> +a=6370000 +b=6370000 +units=m >>>> >>>> +proj=lcc which is usually what wrf is run with. >>>> The tricky part is: >>>> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 >>>> because every WRF run is different (the WRF Preprocessing System >>>> optimizes >>>> the projection for the domain). >>>> and then there is probably no shift so you need(?) +x_0=0 +y_0=0 >>>> >>>> This gives: >>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 >>>> +ellps=sphere >>>> +a=6370000 +b=6370000 +units=m +no_defs >>>> >>>> But, wrf users like to give out lat and long so you need to assign it: >>>> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs >>>> >>>> and then reproject the lat/long to lcc coordinates using this string: >>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 >>>> +ellps=sphere >>>> +a=6370000 +b=6370000 +units=m +no_defs >>>> >>>> One word of caution, make sure you received the correct coordinates. >>>> Some >>>> variables are run cell center while some are run at cell edge. It looks >>>> like from your .nc file it was made by your collaborator so I assume >>>> they >>>> are right. >>>> >>>> That said, another word of caution, I found that the XLAT and XLONG >>>> variables from WRF output aren't very precise. There is a "geogrid" file >>>> from the preprocessing system that has the domain corners, resolution, >>>> nrow >>>> and ncol from which you can make a better grid using the native >>>> projection >>>> system (in my case it was a 4km grid). I suggest you try to get those. >>>> >>>> I hope this helps... I have to run but wanted to save people too much >>>> head >>>> scratching. I can get you running with more help tonight if you need. >>>> Dominik >>>> >>>> >>>> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <agus.cama...@gmail.com> >>>> wrote: >>>> >>>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill >>>>> bring >>>>> more giant crab this time. >>>>> >>>>> The creator of the .nc file also looked at this webpage: >>>>> http://www.pkrc.net/wrf-lambert.html >>>>> It seemed like the right proj4 string might be this one: >>>>> >>>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 >>>>> +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs >>>>> >>>>> However this projection also does not allow me to adequately plot the >>>>> locations on the raster. >>>>> >>>>> Here is the .nc file. it contains several layers. >>>>> >>>>> >>>>> >>>>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0 >>>>> >>>>> >>>>> >>>>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <mdsum...@gmail.com>: >>>>> >>>>> On Tue, 23 Feb 2016 at 20:09 Roger Bivand <roger.biv...@nhh.no> wrote: >>>>>> >>>>>> On Tue, 23 Feb 2016, Alex Mandel wrote: >>>>>>> >>>>>>> I made an attempt at it too. Investigating the original data, I'm >>>>>>>> >>>>>>> not >>>>> >>>>>> sure that the projection information supplied is correct for the >>>>>>>> >>>>>>> data >>>>> >>>>>> linked. When I load up the data in a unprojected space, the >>>>>>>> >>>>>>> coordinates >>>>> >>>>>> don't look at all similar to any Lambert projected data I have, they >>>>>>>> actually look like Lat/Lon in some unprojected coordinate system, >>>>>>>> perhaps a different spheroid than expected. >>>>>>>> >>>>>>> Does anyone have a link to the original data? Is is possible that >>>>>>> >>>>>> this is >>>>> >>>>>> the General Oblique Transformation used by modellers - that is >>>>>>> >>>>>> something >>>>> >>>>>> that feels like longlat but is recentred and oblique? Example at the >>>>>>> >>>>>> very >>>>> >>>>>> end of my GEOSTAT talk last year (slides 81-83): >>>>>>> >>>>>>> http://geostat-course.org/system/files/geostat_talk_150817.pdf >>>>>>> >>>>>>> Roger >>>>>>> >>>>>>> >>>>>>> For what it is worth, the General Oblique Transformation is not the >>>>>> only >>>>>> example - it's very common for modellers to have a mesh that has the >>>>>> "mostly-properties" of a projection, but is not actually describable >>>>>> >>>>> with >>>>> >>>>>> standard transform + affine parameters. The main cases that I've seen >>>>>> >>>>> are >>>>> >>>>>> polar stereographic, equal area or oblique Mercator. Often they really >>>>>> >>>>> are >>>>> >>>>>> simple transforms and you can reconstruct without loss, but it's not >>>>>> usually possible to tell without exploration. It's an interesting >>>>>> dis-connect to see code that builds a mesh with certain properties, >>>>>> then >>>>>> only stores longitudes and latitudes - when it could be done with >>>>>> >>>>> standard >>>>> >>>>>> tools and be stored and used much more efficiently. >>>>>> >>>>>> (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area >>>>>> terminology conflated in this context too. ) >>>>>> >>>>>> I'm also interested to explore the original data. >>>>>> >>>>>> Cheers, Mike. >>>>>> >>>>>> >>>>>> >>>>>> -Alex >>>>>>>> >>>>>>>> On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote: >>>>>>>> >>>>>>>>> Hi >>>>>>>>> >>>>>>>>> I tried to make it work but I had to give up. I wanted to reproject >>>>>>>>> >>>>>>>> the >>>>>> >>>>>>> Lamberth conformal conic coordinates to long-lat but it didn't work. >>>>>>> >>>>>>>> Perhaps someone can see what I did wrong. Here is what I did (data >>>>>>>>> >>>>>>>> in >>>>> >>>>>> R >>>>>> >>>>>>> binary format and figure in png format both attached): >>>>>>> >>>>>>>> library(raster) >>>>>>>>> library(maptools) >>>>>>>>> data(wrld_simpl) >>>>>>>>> >>>>>>>>> r <- raster("raster.grd") >>>>>>>>> projection(r) >>>>>>>>> ## > NA >>>>>>>>> >>>>>>>>> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep = >>>>>>>>> >>>>>>>> ",") >>>>> >>>>>> coordinates(uro) <- ~lon+lat >>>>>>>>> >>>>>>>>> ## Set projections for the 3 data sets >>>>>>>>> >>>>>>>>> ## Lamberth's confocal conic projection with given parameters >>>>>>>>> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0 >>>>>>>>> >>>>>>>> +lat_2=45.0 >>>>> >>>>>> +ellps=WGS84" >>>>>>> >>>>>>>> projection(r) >>>>>>>>> >>>>>>>>> ## Assume that lon, lat are geographical coordinates (degrees >>>>>>>>> >>>>>>>> decimal) >>>>> >>>>>> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84") >>>>>>>>> >>>>>>>>> ## wrld_simpl is in geographical coordinates >>>>>>>>> proj4string(wrld_simpl) >>>>>>>>> >>>>>>>>> ## Make figure in png format >>>>>>>>> ## Of course plotting data with 2 different projections will give >>>>>>>>> ## some distortions >>>>>>>>> pdf("uro.png") >>>>>>>>> >>>>>>>>> plot(r) >>>>>>>>> points(uro) >>>>>>>>> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of >>>>>>>>> >>>>>>>> 'r' >>>>> >>>>>> dev.off() >>>>>>>>> >>>>>>>>> extent(r) >>>>>>>>> ## class : Extent >>>>>>>>> ## xmin : -131.4368 >>>>>>>>> ## xmax : -68.56323 >>>>>>>>> ## ymin : 12.35567 >>>>>>>>> ## ymax : 50.26619 >>>>>>>>> >>>>>>>>> ## Reproject the raster to long-lat >>>>>>>>> ## This doesn't work (collapsed domain) >>>>>>>>> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs >>>>>>>>> >>>>>>>> +ellps=WGS84 +towgs84=0,0,0") >>>>>>> >>>>>>>> ## Because >>>>>>>>> >>>>>>>>>> extent(rp) >>>>>>>>>> >>>>>>>>> ## class : Extent >>>>>>>>> ## xmin : -100.0015 >>>>>>>>> ## xmax : -99.68557 >>>>>>>>> ## ymin : 37.70658 >>>>>>>>> ## ymax : 38.00046 >>>>>>>>> >>>>>>>>> ## Save data in R binary format >>>>>>>>> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData") >>>>>>>>> >>>>>>>>> >>>>>>>>> Yours sincerely / Med venlig hilsen >>>>>>>>> >>>>>>>>> Frede Aakmann Tøgersen >>>>>>>>> Specialist, M.Sc., Ph.D. >>>>>>>>> Plant Performance & Modeling >>>>>>>>> >>>>>>>>> Technology & Service Solutions >>>>>>>>> T +45 9730 5135 >>>>>>>>> M +45 2547 6050 >>>>>>>>> fr...@vestas.com >>>>>>>>> http://www.vestas.com >>>>>>>>> >>>>>>>>> Company reg. name: Vestas Wind Systems A/S >>>>>>>>> This e-mail is subject to our e-mail disclaimer statement. >>>>>>>>> Please refer to www.vestas.com/legal/notice >>>>>>>>> If you have received this e-mail in error please contact the >>>>>>>>> >>>>>>>> sender. >>>>> >>>>>> >>>>>>>>> >>>>>>>>> -----Original Message----- >>>>>>>>> From: R-sig-Geo [mailto:r-sig-geo-boun...@r-project.org] On >>>>>>>>> >>>>>>>> Behalf Of >>>>> >>>>>> Agus Camacho >>>>>>> >>>>>>>> Sent: 22. februar 2016 19:20 >>>>>>>>> To: t...@wildintellect.com >>>>>>>>> Cc: r-sig-geo >>>>>>>>> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a >>>>>>>>> >>>>>>>> reference system implicit in a .nc file >>>>>>> >>>>>>>> Thanks Alex, but the locations still fall in the sea when i plot >>>>>>>>> >>>>>>>> them >>>>> >>>>>> using >>>>>>> >>>>>>>> your recommended Solution. I looked at the sites you proposed and >>>>>>>>> >>>>>>>> they >>>>> >>>>>> have >>>>>>> >>>>>>>> other values for lat_1, lat_0, etc.. >>>>>>>>> >>>>>>>>> 2016-02-22 11:04 GMT-07:00 Alex M <tech_...@wildintellect.com>: >>>>>>>>> >>>>>>>>> On 02/22/2016 09:50 AM, Agus Camacho wrote: >>>>>>>>>> >>>>>>>>>>> Dear all, >>>>>>>>>>> >>>>>>>>>>> Im trying to overlap these points: >>>>>>>>>>> >>>>>>>>>>> >>>>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0 >>>>> >>>>>> and a wrld_simpl object: >>>>>>>>>>> library(maptools) >>>>>>>>>>> data(wrld_simpl) >>>>>>>>>>> >>>>>>>>>>> Over this raster layer >>>>>>>>>>> >>>>>>>>>>> >>>>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0 >>>>> >>>>>> This rastr comes from a .nc file without a reference system. The >>>>>>>>>>> >>>>>>>>>> author >>>>>>> >>>>>>>> of >>>>>>>>>> >>>>>>>>>>> that .nc file gave me the following data about the .nc. >>>>>>>>>>> >>>>>>>>>>> The projection is *Lambert conformal conic* projection >>>>>>>>>>> CEN_LAT = 38.0 >>>>>>>>>>> CEN_LON = -100.0 >>>>>>>>>>> TRUELAT1 = 25. >>>>>>>>>>> TRUELAT2 = 45. >>>>>>>>>>> >>>>>>>>>>> However, despite i have gone through many sites in the internet, >>>>>>>>>>> >>>>>>>>>> i >>>>> >>>>>> cant >>>>>>> >>>>>>>> figure it out: >>>>>>>>>>> >>>>>>>>>>> a) if that is all the data i need to set a reference system for >>>>>>>>>>> >>>>>>>>>> my >>>>> >>>>>> points >>>>>>> >>>>>>>> and the wrld_simp object. >>>>>>>>>>> >>>>>>>>>>> b) how to change a typical CRS object with such data >>>>>>>>>>> >>>>>>>>>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84") >>>>>>>>>>> >>>>>>>>>>> Where do i enter the TRUELAT and CENLAT values? >>>>>>>>>>> Are there any site that explains easily what the fields in the >>>>>>>>>>> >>>>>>>>>> CRS >>>>> >>>>>> mean >>>>>>> >>>>>>>> and >>>>>>>>>> >>>>>>>>>>> how to change them? >>>>>>>>>>> >>>>>>>>>>> Thanks in advance. >>>>>>>>>>> >>>>>>>>>>> https://github.com/OSGeo/proj.4/wiki/GenParms >>>>>>>>>> https://trac.osgeo.org/proj/wiki/GenParms >>>>>>>>>> >>>>>>>>>> I believe: >>>>>>>>>> +lat_0 = CEN_LAT Latitude of origin >>>>>>>>>> +lat_1 = TRUELAT1 Latitude of first standard parallel >>>>>>>>>> +lat_2 = TRUELAT2 Latitude of second standard parallel >>>>>>>>>> +lon_0 = CEN_LON Central meridian >>>>>>>>>> >>>>>>>>>> proj strings are defined by the proj4 libary. It's website listed >>>>>>>>>> >>>>>>>>> above >>>>>> >>>>>>> and the associated mailing lists or gis stackexchange would be the >>>>>>>>>> places to get help on it. >>>>>>>>>> https://lists.osgeo.org/mailman/listinfo/metacrs >>>>>>>>>> >>>>>>>>>> It often helps to browse similar projections on >>>>>>>>>> http://spatialreference.org/ >>>>>>>>>> http://epsg.io/ >>>>>>>>>> >>>>>>>>>> Enjoy, >>>>>>>>>> Alex >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>> _______________________________________________ >>>>>>>> R-sig-Geo mailing list >>>>>>>> R-sig-Geo@r-project.org >>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>>>>>> >>>>>>> -- >>>>>>> Roger Bivand >>>>>>> Department of Economics, Norwegian School of Economics, >>>>>>> Helleveien 30, N-5045 Bergen, Norway. >>>>>>> voice: +47 55 95 93 55; fax +47 55 95 91 00 >>>>>>> e-mail: roger.biv...@nhh.no >>>>>>> http://orcid.org/0000-0003-2392-6140 >>>>>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en >>>>>>> http://depsy.org/person/434412 >>>>>>> _______________________________________________ >>>>>>> R-sig-Geo mailing list >>>>>>> R-sig-Geo@r-project.org >>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>>>>> >>>>>> -- >>>>>> Dr. Michael Sumner >>>>>> Software and Database Engineer >>>>>> Australian Antarctic Division >>>>>> 203 Channel Highway >>>>>> Kingston Tasmania 7050 Australia >>>>>> >>>>>> [[alternative HTML version deleted]] >>>>>> >>>>>> _______________________________________________ >>>>>> R-sig-Geo mailing list >>>>>> R-sig-Geo@r-project.org >>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>>>> >>>>>> >>>>> >>>>> -- >>>>> Agustín Camacho Guerrero. >>>>> Doutor em Zoologia. >>>>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de >>>>> Biociências, USP. >>>>> Rua do Matão, trav. 14, nº 321, Cidade Universitária, >>>>> São Paulo - SP, CEP: 05508-090, Brasil. >>>>> >>>>> [[alternative HTML version deleted]] >>>>> >>>>> _______________________________________________ >>>>> R-sig-Geo mailing list >>>>> R-sig-Geo@r-project.org >>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>>> >>>>> >>>> >>> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> >> > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Agustín Camacho Guerrero. Doutor em Zoologia. Laboratório de Herpetologia, Departamento de Zoologia, Instituto de Biociências, USP. Rua do Matão, trav. 14, nº 321, Cidade Universitária, São Paulo - SP, CEP: 05508-090, Brasil. [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo