Monday 31 October 2016

How to find polygons inside polygon with PostGIS 2.0


I have a relation representing all the counties in my country and the other one containing all the districts - e. g. each and every county is composed of several districts. How do I select them with spatial query? I tried with


select nazorp from kraje as k, orp_wgs as o WHERE ST_Within(o.geom,k.geom) AND
k.nazev = 'Liberecký'


but had wrong results returned (well, they are probably not wrong, they are just not what I expected them to be). ST_Within only returns districts that don't share a boundary with the county, but I need to get all the districts within the county. Is that possible? I haven't found any built-in function suitable for my needs yet.




arcgis desktop - Changing raster cells values inside of a polygon to the closest cells just outside the polygon!


I have been wrestling with this for a while and couldn't find out how I can do this. I want each cell of a raster inside of a polygon boundary to have the value of the closest cell just outside of the polygon.



I think the closest answer was the one at this question, but I couldn't apply it!


My plan was:



  1. Use "Polygon to Raster" to change the polygon to a raster. The cells inside the polygon become 1 and the points outside of the polygon are nodata.

  2. Use Con(IsNull(Raster),0,NULL) and assign it to Raster2

  3. Add Raster1 and Raster2 to get nodata inside the polygon

  4. It should be easy from here with nodata inside the polygon. We just use "Euclidean Allocation" and it does the job.


I'm stuck at step 2 and get an error:




[SELECT * FROM VAT_Extent_RAS WHERE IsNull( "Value" )] Failed to execute (Con).



Do you have any thoughts? Is there an easier way to do this?


Update:


I'll try to summarize my progress here. I want to change all raster cells inside the polygon to "nodata" while the raster cells outside of the polygon remain the same. I used polygon to raster function first. This gives me a raster with 0 cell inside the polygon and "nodata" outside of the polygon. What do I do next?




arcgis server - Toggle Toolbar Button and Activate Identify


I already asked how to toggle an identify button in a previous question, and was given the correct answer to the question, but when I tried to implement it into a toolbar that would disable the identify button when another button was pressed it didn't work. I guess I should have asked the question a little better. I am hoping that soneone will be able to help put the two pieces of code together. An example can be found at the following link and my code is below.


http://geoville.org/viewers/PopUp/


    







Identify with Popup






















Thanks to those that can help and thanks to those that helped with the previous question.


-Mike



Answer



Ok, you are needing to add some additional logic to know what other buttons are activated and toggle on and off as needed. Something like this:


    function updateTool(theDijit) {
var zoominDij, zoomoutDij, panDij, idDij, itemArray;
zoominDij = dijit.byId("zoomin");
zoomoutDij = dijit.byId("zoomout");

panDij = dijit.byId("pan");
idDij = dijit.byId("identify");
itemArray = [zoominDij, zoomoutDij, panDij, idDij];
if (theDijit.attr("id") != "identify")
identifyDisconnect();
for (var i = 0; i < itemArray.length; i++) {
if (theDijit.attr("id") == itemArray[i].attr("id")) {
theDijit.attr("checked", true);
} else {
itemArray[i].attr("checked", false);

itemArray[i].focus();
}
}
}

You will need a little work on your itemArray to have the right button IDs as well as inside where you set the trigger to disable the other functions.


Coupled with a call like this on your toolbar buttons:


onclick="identifyConnect();updateTool(this);"

Now this additional call will set a value of what button was clicked, will help you disconnect the function from one and allow you to set active to the other. This is what I use on my productions sites like here Washington State Sales Tax.



geoserver - WFS Request Layer (GML) with Filter


Is there a solution to request the GML layer from Geoserver with filter like using bbox:


http://localhost:8080/geoserver/tiger/ows?service=WFS&version=1.0.0
&request=GetFeature&typeName=tiger:poi&maxFeatures=50&
bbox=-74.0104611,40.70758763,-74.00153046439813,40.719885123828675

The result will be specific filtering the bounding box (lat/lon), based on bbox parameter that user entered.


The Result:





unknown






-74.00857344,40.71194565


lox
pics/22037884-Ti.jpg
pics/22037884-L.jpg





The problem now, I want to filter the data based on the lox, already tried with entered the tiger:NAME=lox as a parameter when request, the result just same like when enter no parameter. Is there a solution for this?



Answer



Based on Filter Encoding with Spatial Filter within WFS-Request and Geoserver Filter. I found that there is a parameter named filter and the value can enter like this:


http://localhost:8080/geoserver/tiger/ows?service=WFS&version=1.0.0
&request=GetFeature&typeName=tiger:poi&maxFeatures=2
&filter=NAME
lox

Seeking open source software similiar to Esri City Engine?




I am looking for open source software similar to ESRI City Engine and currently going through open source software directory.


About ESRI City Engine (Quick introduction)


We want to make few historical places in 3D format & I heard that in some cases its difficult to integrate with other software (i.e. Esri City engine to Google earth and vice verse etc)



Most of the spatial data is in PostgreSQL/PostGIS . Hence looking for alternative so any open source software included 3 D editing facility will be great.



Answer



You could look at Blender which is the defacto standard open source 3D modeling program with a huge community and tons of tutorials out on the web.


Here is a tuturial about creating a cityscape in Blender.


modelling - Moving window regression (MWR) in R?


I'm looking for some package in R to estimate the value of a variable by moving window regression. I have read in Local Models for Spatial Analysis, which is a good method to predict values ​​at unsampled locations. A very similar case is mentioned in the book, PRISM.


I've been reading various R packages, but have not found anything about MWR, GWR only.


Is there another way of naming or can GWR be adapted to function as MWR?


Text from book -> Local Models for Spatial Analysis




Combine Leaflet with Geoserver WMS map



I'm not very skilled within the GIS world, however, I would like to know if it's possible to use this map:


http://wfs-kbhkort.kk.dk/k101/wms?service=WMS&version=1.1.0&request=GetMap&layers=k101:theme-startkort&styles=&bbox=12.451709828167994,55.63150896601625,12.64714535243575,55.72503797474024&width=512&height=495&srs=EPSG:4326&format=application/openlayers


within leaflet?


I've looked at the API for leaflet and can see that it is possible to use custom maps/tiles with the tileLayer option. However, I cannot get it to work.


Can someone point me in the right direction?



Answer



Take a look at L.tilelayer.wms in the API.


In your case you should get something like this:


var map = L.map('map').setView([55.67, 12.60], 11);


var mywms = L.tileLayer.wms("http://wfs-kbhkort.kk.dk/k101/wms", {
layers: 'k101:theme-startkort',
format: 'image/png',
transparent: true,
version: '1.1.0',
attribution: "myattribution"
});
mywms.addTo(map);

shapefile - QGIS: Release file lock on file used in a processing algorithm at the end of script


I am using runalg (qgis:joinattributesbylocation) on two QgsVectorLayer objects. At the end of my script I'm trying to delete all the files created during the script. I can delete everything except the .shp and .dbf files from the two input vector layers for the process algorithm above. I remove all layers from the QgsMapRegistry before I try to delete the files. I've tried deleting all variables pointing to the layers, still doesn't work. The file lock remains until I exit QGIS. I tried the workaround posted here and it doesn't work for me, re-adding and removing the layers from the registry doesn't work. I also created new layers with the same files, added those to the registry, and then remove those, still no luck. The only way I can delete the files is to close QGIS.





Sunday 30 October 2016

Postgis Self-intersection on geometryraster comparison


I am rather new to postgis and try to get a (weighted) average height of an area around a point. I am running Postgresql 9.5.5 with Postgis 2.2.1 on Ubuntu. I have imported ASTER GDEM elevation information as a raster and built the following functions, following this tutorial:


WITH statics AS (
SELECT (ST_DISTANCE(CAST (ST_Project(CAST(ST_SETSRID($1, 4326) AS geography), $2, radians(90.0)) AS geometry), CAST(ST_SETSRID($1, 4326) AS geometry))) AS exp_fact
-- calculate degrees to expand the point with, $2 is given in meters
)
SELECT SUM(sub.height * sub.weight) / SUM(sub.weight)

FROM (
SELECT height, ((1 - (dist / $2))) as weight
FROM get_elev_points_within_geom(ST_MAKEVALID(ST_BUFFER(ST_SETSRID($1, 4326), (SELECT exp_fact from statics))), $1)
WHERE dist <= $2
) as sub

where get_elev_points_within_geom(geo geometry, center geometry) is:


SELECT ST_Y(ST_CENTROID((a.neighborhood).geom)),
ST_X(ST_CENTROID((a.neighborhood).geom)),
(a.neighborhood).val,

ST_DISTANCE(CAST (ST_CENTROID((a.neighborhood).geom) AS geography), CAST ($2 AS geography)) as dist_m
FROM(
SELECT ST_INTERSECTION(elevation.rast, $1) AS neighborhood
FROM elevation
WHERE ST_INTERSECTS(elevation.rast, $1)
) AS a
ORDER BY dist_m
LIMIT 1000;

Running e.g. this query:



SELECT * FROM get_elev_average_weighted_by_dist_at(ST_MAKEPOINT(8.4305360482003504, 49.013877361340597), 50);

I get:


NOTICE:  Ring Self-intersection at or near point 8.4315277777777791 49.016527777777782
CONTEXT: PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
SQL function "get_elev_points_within_geom" statement 1
SQL function "get_elev_average_weighted_by_dist_at" statement 1
Total query runtime: 145 msec
1 row retrieved.


I have traced the Problem down to to:


SELECT (ST_INTERSECTION(elevation.rast, ST_SETSRID(ST_GeomFromText('POLYGON((8.4298 49.0131,8.4298 49.0145,8.4312 49.0145,8.4312 49.0131,8.4298 49.0131))'),4326))).geom
FROM elevation
WHERE ST_INTERSECTS(elevation.rast, ST_SETSRID(ST_GeomFromText('POLYGON((8.4298 49.0131,8.4298 49.0145,8.4312 49.0145,8.4312 49.0131,8.4298 49.0131))'),4326))
LIMIT 14;

which will return with the same error message (but different trace, obviously). Interestingly, if I reduce the limit to 13, I will not get the Message.


Also this will run for about 150ms per Query on my machine. This is quite annoying if I query for lots of points.


How can I improve on efficiency?





Labeling line features using PyQGIS?


I used the following code to label the line features(roads) in a layer.It works in the console, but not by python code(application) It works for the points too.


What did I do wrong?


 def labell(self):
layer = self.iface.mapCanvas().currentLayer()
palyrr = QgsPalLayerSettings()

palyrr.readFromLayer(layer)
palyrr.enabled = True
palyrr.fieldName = 'name'
palyrr.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'08','')
palyrr.writeToLayer(layer)
layer.commitChanges()
self.iface.mapCanvas().refresh()

Answer



I don't know what could be the difference between console and plugin, I think it should to work either way. Also you try to replace layer.commitChanges() (which makes not sense there) with layer.triggerRepaint() and look if that helps (note it is not necessary to call mapCanvas.refresh()). You can also use custom properties for the layer to enable labeling, a sample code is the following:


layer = iface.activeLayer()

layer.setCustomProperty("labeling/fieldName", "name" )
layer.setCustomProperty("labeling/placement", QgsPalLayerSettings.Line)
layer.setCustomProperty("labeling/fontSize","8" )
layer.setCustomProperty("labeling/enabled","true" )
layer.triggerRepaint()

EDIT


for linear features you have to set the placement property as well. I changed the above code. Your first posted code should also work by adding the following line:


palyrr.placement = QgsPalLayerSettings.Line


You can use other placement attributes please have a look at this link


Saturday 29 October 2016

tomcat - Missing 'Spec' Parameter using print.pdf on GeoServer Print Module


I have two configurations of Tomcat with the same version of GeoServer and the print module installed as below.



  • The old configuration is Tomcat 8.0.44/GeoServer 2.11.1

  • The new configuration is Tomcat 9.0.10/GeoServer 2.11.1


Using the old configuration I can POST to localhost:1080/geoserver/pdf/print.pdf in SoapUI with the spec file generated from the Print Module demo in the body and this returns a pdf. An application has been written in C# based around this method of generating a pdf.


Using the new configuration (as Tomcat 8.0.x is now end of life, and also a combination of Tomcat 8.5.32/GeoServer 2.11.1) the same POST returns an error message:



generating PDF:
Missing 'spec' parameter
]]>


If the POST is changed to a GET with ?spec={SPEC} then this works but for our actual application requirements a POST is the only way to ensure that all of the vector data gets sent to the print module.


Does anyone have experience of using Tomcat 8.5 or above with GeoServer's Print Module using a POST to print.pdf?


The default SPEC file that has been used in SOAPUI testing is


{
"units":"degrees",
"srs":"EPSG:4326",

"layout":"A4 portrait",
"dpi":75,
"layers":[
{
"baseURL":"http://localhost:1080/geoserver/wms",
"opacity":1,
"singleTile":false,
"type":"WMS",
"layers":["topp:states"],
"format":"image/jpeg",

"styles":[""],
"customParams":{}
},
{
"type":"Vector",
"styles":
{
"1":
{
"externalGraphic":"http://geoserver.org/img/geoserver-logo.png",

"strokeColor":"red",
"fillColor":"red",
"fillOpacity":0.7,
"strokeWidth":2,
"pointRadius":12
}
},
"styleProperty":"_gx_style",
"geoJson":{
"type":"FeatureCollection",

"features":[{
"type":"Feature",
"id":"OpenLayers.Feature.Vector_52",
"properties":{"_gx_style":1},
"geometry":{
"type":"Polygon",
"coordinates":[[[-97,39],[-98,40],[-96,41],[-97,39]]]
}
},
{

"type":"Feature",
"id":"OpenLayers.Feature.Vector_61",
"properties":{"_gx_style":1},
"geometry":{
"type":"LineString",
"coordinates":[[-97,40],[-98,39],[-99,38]]
}
},
{
"type":"Feature",

"id":"OpenLayers.Feature.Vector_64",
"properties":{"_gx_style":1},
"geometry":{
"type":"Point",
"coordinates":[-98,38]
}
}
]},
"name":"vector",
"opacity":1

}],
"pages":[{
"center":[-98.000000000002,40],
"scale":4000000,
"rotation":"0",
"mapTitle":"A custom title",
"comment":"A custom comment"
}]
}


coordinate system - Influence of the scale factor on the projection


I'm currently reading documentation about map projections to understand the source code of the Proj4 project.


The scale factor is named in a variety of sources I read. This sources explained its definition and its value for some projections.


In the source code of Proj4, for the mercator projection (sphere and ellipse case), the scale factor influences the coordinates on the projection :


//P->k0 is the scale factor
xy.x = P->k0 * lp.lam;
xy.y = - P->k0 * log(pj_tsfn(lp.phi, sin(lp.phi), P->e));

Why and how I should use the scale factor during the computation of the projection ? Is there any valuable resources on the web ?



This question in asked in the sense of projection computation. I can find the formula for the inverse and forward projection as well as the scale factor in a various of resources, but no one explains how I should use both in a algorithm. You have the definition of the projection and the definition of the scale factor, but it's not clearly written that I should multiply or divide the result by the scale factor.


Is it a general rule : if I find the formula for any projection with the related scale factor, should I, in all cases, always divide or multiply the results by the scale factor ?



Answer



Unfortunately, the term "scale factor" is ambiguous. In cartography, maps and projections, the concept and application of scale is of fundamental importance. By definition, it is a factor – meaning it is something to multiply or divide – so whether "scale" or "factor" is the adjective the particular word pairing has no obvious meaning, except in a particular context:



The map or globe scale is a ratio of distance on map or globe to corresponding distance on the ground or reality. Either it has no units or its units reflect the map and ground units – miles per inch, km per cm, etc. It is variously called scale, map scale, principal scale, representational fraction, or nominal scale. I like that last one, nominal scale, because most maps have a single statement of its scale. Sadly, it is sometimes also called scale factor.



All map projections distort linear scale, all over the map. This distortion is almost always termed scale factor (and sometimes "projection scale factor" or "point scale factor"). At any point on the map, it is the ratio of the "true" (undistorted) scale and the "nominal" (distorted) scale. In other words, it is the ratio of the true ground distance to the implied distance on the map.



When computing a map projection, that is



(X, Y) ← projection (λ, φ)


you need to have some constant which depends on the size of your region of interest, the size of your map, and the map units involved. That constant is our friend the nominal scale. Since you don't provide the full code, and it's a little cryptic, I cannot say for certain, but I suspect that is what is meant by "scale factor" in your particular problem.


According to Mathematics_of_the_Mercator_projection on wiki, the spherical case, which uses radius, R, as a substitute for scale:


X = R λ


Y = R ln [tan (π/4 + φ/2)]


(That looks similar to your code.)


How is radius a substitute for scale? Simple. It is the constant which determines the size of the map: a larger globe yields a larger map. If R is the Earth's radius, then your map scale will be one-to-one. If R is the radius of your globe in, say, mm, inches, or pixels, then X,Y will be in those same units and the map's nominal scale, NS, will be the ratio of your globe's radius to the Earth's radius:


NS = globe radius / Earth radius


To get ground distances from a measurement on a projected map, including any distortions:


ground distance = map distance / NS



To remove distortion, see below.



To properly correct distortion at any point on a projected map, you ought to be able to calculate the distortion, SF (scale factor):


SF ← distortion (λ, φ)


In this case, SF has to be calculated, or provided in a look-up table, wherever it is needed. Does your code calculate "scale factor" as a function of "lam" and "phi"? I doubt it.


According to Mathematics_of_the_Mercator_projection on wiki, which uses K for scale factor:


Spherical case: K = sec φ


Ellipsoidal case: K = sec φ sqrt(1 - e2 sin2 φ)


where e2 is about 0.006 for all reference ellipsoids.


To correct for any projection distortion, i.e., to convert a projected map distance to get true (undistorted) globe distance, always divide by the scale factor, SF:



ground distance = map distance / SF


That might look familiar.



Yes, they're used in exactly the same way. However, whenever the Earth is really projected onto a map that is actually measured in terms of map units (mm, cm, inches, pixels, etc.), then you need to apply both



  • a global nominal scale to get the correct globe/earth magnitude and units, and

  • a local scale factor to correct for projection distortion at any particular Earth position.


If you are not making measurement on an actual map, and you are just using coordinates that are in the same units as your Earth radius, R, then your nominal scale is trivial (NS = 1) and you only need to use the scale factor.


javascript - leaflet count by featuretypes in featuregroup


Is there any way in leaflet to get the number of rectangles in a feature group? I know this code that work: drawnItems.getLayers().length; but it counts all the objects within the drawnItems featuregroup. I only need the a specific type of features, for example rectangle



Answer




Assuming that you are using instances of L.Rectangle, you can use functional programming like so:


drawnItems.getLayers().filter( function(l) {return l instanceof L.Rectangle} ).length;

The filter method of Array will return another array, containing only the elements which make the callback function return true; in this case, it will return an array containing only those layers which are instances of L.Rectangle.


pyqgis - How to rotate a SVG marker via a data-defined override to the value of one field in a layer


I have a pyqgis stand alone app that display linear features. I put an arrow svg symbol (QgsSvgMarkerSymbolLayerV2) on the line, now I would like to rotate it based on the 'value' field of each line feature in the layer. I tried some variations to the suggested approach below, but they are not working:


sl = QgsSvgMarkerSymbolLayerV2('arrow.svg')
lDD = QgsDataDefined(True, True, 'CASE WHEN "value" < 0 THEN 180.0 ELSE 0.0 END CASE', 'value')
sl.setDataDefinedProperty('angle', lDD)



How do I plot points as graduated/proportional circles in R?


I'm back again with more R mapping questions! This one is concerning symbolizing points by a certain attribute with graduated/proportional circles.


DATA: My full code can be found here if you're interested. My CSV of interest is here. Ethiopia shapefile can be downloaded here


I've plotted points from the CSV over a shapefile. For the sake of brevity, here's the heavily abridged version of my code:


library(raster)
library(rgdal)

#set your working directory


#read in eth shapefile
eth <- readOGR(dsn = "D:/Mapping-R/Ethiopia", layer = "ETH_adm0")

#read in total branch location csv
branches <- read.csv("Branches_Africa.csv", header = TRUE)

#transform branches data frame to spdf for mapping
coordinates(branches) <- ~Lon + Lat

#plot shapefile with point overlay

plot(eth)
points(branches$Lon, branches$Lat, col = "blue", pch = 16, cex = .5)

As you can see in the branches CSV, there are several different attributes. In my case, I would like to symbolize the "share" attribute with graduated/proportional circles, as seen here: enter image description here


Ideally I would also like to change the colors of each point as well to match the QGIS symbology as closely as possible.


If you look at my full code, you'll notice that I've not used ggplot to map, so I would prefer the solution not come from that package (or would ggplot make things easier once I've transformed all the data to fit the package?). I tried using spplot as per this post with:


spplot(branches, "share", col.regions = brewer.pal(9, "Reds"), cuts = 8, scales = list(draw = T))

but the code gave me the following error: Error in fill.call.groups(args.xyplot, z = z, edge.col = edge.col, colorkey = colorkey, : number of colors smaller than number of factor levels


I also tried to work with mapCircles from the rCarto package, but I had no luck there either.



SO, how can I change the symbology of branches$share so that graduated/proportional (and preferably differently colored) circles are displayed instead of just one size/color?



Answer



I know you prefer the solution not to use the ggplot package, but I thought it might be helpful to add a quick example, just for comparison.


library(rgdal)
library(ggmap)
library(scales)

#data
eth <- readOGR(".", "ETH_adm0")
branches <- read.csv("Branches_Africa.csv", header = T)

branches <- branches[branches$CO == "ET", ] # select Ethiopia subset
branches$share <- as.numeric(as.character(branches$share)) # convert to numbers

# fortify Ethiopia boundary for ggplot
eth_df <- fortify(eth)
# get a basemap
eth_basemap <- get_map(coordinates(eth), zoom = 5)

# plot
ggmap(eth_basemap) +

geom_polygon(data=eth_df, aes(x=long, y=lat, group=group), fill="red", alpha=0.1) +
geom_point(data=branches, aes(x=Lon, y=Lat, size=share, fill=share), shape=21, alpha=0.8) +
scale_size_continuous(range = c(2, 9), breaks=pretty_breaks(7)) +
scale_fill_distiller(breaks = pretty_breaks(7))

enter image description here


Join spatial point data to polygons in R


I am trying to perform a spatial join between point data and polygon data.


I have data that indicate the spatial coordinates of an event in my csv file A and have another file, shapefile B, that contains the boundaries of an area as polygons.


head(A)
month longitude latitude lsoa_code crime_type
1 2014-09 -1.550626 53.59740 E01007359 Anti-social behaviour
2 2014-09 -1.550626 53.59740 E01007359 Public order
3 2014-09 -1.865236 53.93678 E01010646 Anti-social behaviour


head(B@data)
code name altname
0 E05004934 Longfield, New Barn and Southfleet
1 E05000448 Lewisham Central
2 E05003149 Hawcoat

I want to join the crime data A to my shapefile B to map the crime events that happen in my area A. Unfortunately I cannot perform an attribute join based code as the code in A refers to different units than the code in B.


I've read a number of tutorials and posts but could not find an answer. I tried:


joined = over(A, B)


and overlay, but did not accomplish what I wanted.


Is there a way to do this join directly or would an intermediate transformation from A to another format be needed?


Conceptually I want to select those points of A that fall into the code areas of B (similar to "join based on spatial location in ArcGIS").


Did someone have this issue and solved it?



Answer



The point.in.poly function in the spatialEco package returns a SpatialPointsDataFrame object of the points that intersect an sp polygon object and optionally adds the polygon attributes.


First lets add the require packages and create some example data.


require(spatialEco)
require(sp)

data(meuse)
coordinates(meuse) = ~x+y
sr1=Polygons(list(Polygon(cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409,
180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676,
332618, 332413, 332349)))),'1')
sr2=Polygons(list(Polygon(cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437,
179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
331133, 331623, 332152, 332357, 332373)))),'2')
sr3=Polygons(list(Polygon(cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875,
179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),

c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
329783, 329665, 329720, 329933, 330478, 331062, 331086)))),'3')
sr4=Polygons(list(Polygon(cbind(c(180304, 180403,179632,179420,180304),
c(332791, 333204, 333635, 333058, 332791)))),'4')
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data.frame(row.names=c('1','2','3','4'), PIDS=1:4, y=runif(4)))

Now, lets take a quick look at the data and plot it.


head(srdf@data)  # polygons
head(meuse@data) # points

plot(srdf)
points(meuse, pch=20)

Finally, we can intersect the points with the polygons. The results will be a SpatialPointsDataFrame object with, in this case, two extra attributes (PIDS, y) that were contained in the srdf polygon data.


  pts.poly <- point.in.poly(meuse, srdf)
head(pts.poly@data)

If there is not a unique identification column in the polygon data you could easily add one.


srdf@data$poly.ids <- 1:nrow(srdf) 


Once we have the points and polygons intersected, we can aggregate the points using the unique polygon ID's that were an attribute in the polygon data.


# Number of points in each polygon
tapply(pts.poly@data$lead, pts.poly@data$PIDS, FUN=length)

# Mean lead in each polygon
tapply(pts.poly@data$lead, pts.poly@data$PIDS, FUN=mean)

Editing date in *.ssf files?



So my trimble TSC1 Pro XRS (bought in 2000) is getting old I guess and it change byitself the date. I collected data this week, but it indicate from 1998.


I view this info because I can't apply diferential correction. Indeed, the post process file from 1998 aren't avaible anymore.


So I want to change the change the date of the record in *.ssf files.


I tried two different ways : - The office pathinder (4.0) SSF record editor. It doesn't allow to edit the date.



  • I copy and paste all the record from the *.ssf files in notepad ++, then change the date, then change the files extention in *.ssf. But when I open this files in office pathinder, it doesn't open it correctly.


I'm still looking for any option.



Answer



I finally found a solution :




  • Open *.ssf files in office pathinder

  • Menu tools> other > check *.ssf files

  • Select the right files

  • Options : check modify GPS week number manually

  • Fill wiht the right GPS week number, that you going to find here : https://www.labsat.co.uk/index.php/en/gps-time-calculator

  • Then you can apply the right post-traitement.


Otherwise my date problem seems to come from my TSC1.


arcgis desktop - How much RAM can ArcMap use?


I am exporting 36"x48" label intensive maps as PDF. For urban areas, ArcMap 10 hangs or kicks back an error saying it could not complete the operation. I have a quad processor with 4 GB ram. I also increased cache size to 100GB and have a large pagefile size, etc.


Will increasing physical RAM to 6GB help, or is this a software limit to how much labeling it can handle?



Answer



http://support.esri.com/en/knowledgebase/techarticles/detail/38343 seems to cover it. Basically, it's still a 32 bit application, so it depends on the process whether you can squeeze more than 2 gb of usage out. However, more RAM in the system means less competition for that 4 GB by all the other processes.


Friday 28 October 2016

postgis - Can Geoserver return the raster value of a lat/lon point


I have a GeoTIFF published by Geoserver 2.2 and being used by Leaflet via WMS. How can I get Geoserver to return the raster value at a lat/long point?


I'm guessing it has to do with using WCS, but no examples seem to exist for doing this!


This will be similar to querying PostGIS with ST_Value(raster, ST_SetSRID(ST_Point(lon,lat),4326)). I dont really want to import the same raster into PostGIS just to find the raster value at a point. Or is this the recommended solution?




qgis - Why is my new Shapefile layer not showing up?


I am having trouble with a new shape layer not showing up on my map. The attribute table is showing up fine and all of the points are in there, however visually, there are no dots showing on the map. I have made sure that the new shape layer was at the very top of the stack. I tried placing the new shape layer on a blank new project and it worked perfectly as it should. Then I re-opened the original problem project, deleted all of the existing layers, saved it as temp with no data in it at all, and then closed and re-opened the project and re-loaded the problem shape layer into the new temp project. It did not show up again.


I have gone back to my original data and re-created the csv that I initially used to create the shape file (did this a couple of times), in case there was some spurious content in there, however, this had no effect, and the problem remains.


Using QGIS 1.8.0 and Windows 7 64-bit. However I installed QGIS onto a Linux machine and the identical problem occurred.


Hope someone can help, or has an idea what it might be...




Answer



if you got error in both OS and you can see the information on attribute table, maybe there is a problem with the Coordinate Reference System (CRS).



  • First, when you created the shapefile through "Add Delimited Text Layer" tool, did you assign the CRS?;

  • Second, when you create the project, are you sure that all layers are in the same CRS; If do not, did you activated the option "Enable on the fly CRS tranformation" on the Project Properties window?


CRS tab


pyqgis - Programmatically exporting composition as image from QGIS?


QGIS Print Composer has tool "Save As Image" and I found appropriate function in QGIS source code exportCompositionAsImage. Is it possible to get access to existing Print Composer from Python console and perform this operation - export composition as image with composer settings.


UPDATE #1: Thanks to @ndawson. I'm using printPageAsRaster function for printing different combinations of layers created with layerCombinations plugin:



from layerCombinations import classFactory
lc = classFactory(iface)
lc_manager = lc.manager
combinations = lc_manager.combinationsList
composition = iface.activeComposers()[0].composition()
composition_items = composition.items()
for cmb in combinations:
for item in composition_items:
if item.type() == 65641: #this is the type of ComposerMap
lc_manager.applyCombinationToMap(cmb, item)

image = composition.printPageAsRaster(0)
image.save('/home/dr/' + cmb +'.png','png')

Answer



To get a list of existing composer windows:


composers = iface.activeComposers()

Then, you need to get a reference to a specific composition:


c = composers[0].composition()

Render the page as an image (0 refers to the page number):



image = c.printPageAsRaster(0)

Save the image to disk (first argument is the filename, second is the format):


image.save('output.png','png')

So, the entire code should look something like this:


c = iface.activeComposers()[0].composition()
image = c.printPageAsRaster(0)
image.save('output.png','png')

arcgis server - arcpy.PackageWorkspace doesn't keep reference to SDE enterprise geodatabase


I am publishing a geoprocessing service on ArcGIS Server v10.3.0 that updates data on an ArcSDE Enterprise Geodatabase.


I have the following at the beginning of my script:


arcpy.env.workspace = "Database Connections/TestDatabase@hq-gissql01.sde"


As discussed at Connecting Geoprocessing Service to ArcSDE?, when I publish the script as a geoprocessing service, that line gets changed to:


arcpy.env.workspace = arcpy.env.packageWorkspace

This is all well and good, but arcpy.env.packageWorkspace is simply the v101 folder for this geoprocessing task. It is not the actual .SDE connection file. I can confirm that the .SDE connection file does get moved to the v101 folder, but arcpy.env.packageWorkspace is simply a file path to the v101 folder--not that file.


Because of this, when I call arcpy.ListFeatureClasses() or arcpy.ListDatasets(), an empty list gets returned. From what I can tell this is because the workspace is set to a folder on the server's file system, and not an .SDE file.


How are .py files being re-written when they are published as geoprocessing services? What is the proper syntax to use when setting the workspace so that the published version of the script references an .SDE file and not simply a file folder?



Answer



Try this...change arcpy.env.workspace = "Database Connections/TestDatabase@hq-gissql01.sde" to


import os myfolder = 'c:\mysdefiles' sdePath = os.path.join(myfolder, 'TestDatabase@hq-gissql01.sde') arcpy.env.workspace = sdePath


Then, you obviously need to copy the sde file from the database connections node to the folder you've referenced in the myfolder variable. Also register that folder with the servers data store.



WinError 32 error when trying to delete shapefiles with QGIS python


I've looked at some answers on here about deleting shapefiles with QGIS but I still seem to be running into errors. I am trying to delete shapefiles that have no attribute data in them (so 0 rows). My main issue is that I get:


PermissionError: [WinError 32] The process cannot access the file because it is being used by another process:


I have tried using QgsVectorFileWriter.deleteShapeFile(f) but that leaves the .shp and .dbf behind.


I have also tried using os.remove for each individual file extension and that gives me the Win Error 32.


It seems like my files are still being used in QGIS. Does anyone know how to work around it?



Here is my full script:


clipped_soilpoly = get_data(clipped_folder, ".shp") # makes a list of all the .shps
for f in clipped_soilpoly:
shapelayer = QgsVectorLayer(f,"clipped_poly")
rowcount = shapelayer.featureCount()`

if rowcount < 1:

print ("deleting " + f + " - does not intersect")


#QgsVectorFileWriter.deleteShapeFile(f) <-- only deleting some of the file extensions

split_path = os.path.splitext(f)[0]
del f

#delete .shp
shp = split_path + ".shp"
os.remove(shp)

#delete .dbf

dbf = split_path + ".dbf"
os.remove(dbf)

#delete .prj
prj = split_path + ".prj"
os.remove(prj)

#delete .qpj
qpj = split_path + ".qpj"
os.remove(qpj)


#delete .shx
shx = split_path + ".shx"
os.remove(shx)


else:
print ("keeping " + f)

Answer



On Windows you must stop using and close the file before you can delete it. So QgsVectorFileWriter.deleteShapeFile(f) will work, once you have let go of the file which is still being used by shapelayer.



The QgsVectorLayer is a wrapper around an OGR C++ call so the easiest way to dispose of it is to set it to None.


clipped_soilpoly = get_data(clipped_folder, ".shp") # makes a list of all the .shps
for f in clipped_soilpoly:
shapelayer = QgsVectorLayer(f,"clipped_poly")
rowcount = shapelayer.featureCount()

if rowcount < 1:
shapelayer = None
QgsVectorFileWriter.deleteShapeFile(f)


should work.


raster - Convert x, y position in georeferenced image (with world file) to longitude, latitude?


What steps do I need to do in order to calculate longitude and latitude for a specific pixel in a georeferenced image? For example this image http://www.shadedreliefarchive.com/medit_Mideast_africa_copy.html contains information, .tfw file:


1627.32174969982000 
0.00000000000000
0.00000000000000
-1627.32174969982000

-4957506.34276705000000
6096922.76188281000000

and .prj file:


PROJCS["Lambert Azimuthal Equal-Area",
GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984", SPHEROID["WGS_1984", 6378137, 298.257223563]],
PRIMEM["Greenwich", 0],
UNIT["Degree", 0.017453292519943295]],
PROJECTION["Lambert_Azimuthal_Equal_Area"],

PARAMETER["false_easting", 0],
PARAMETER["false_northing", 0],
PARAMETER["latitude_of_origin", 0],
PARAMETER["central_meridian", 20],
PARAMETER["xy_plane_rotation", 0],
UNIT["Meter", 1]]

I've created this snippet of javascript to compute new x1 and y1 values for specific pixels in the image:


var A = 1627.32174969982000,
D = 0.00000000000000,

B = 0.00000000000000,
E = -1627.32174969982000,
C = -4957506.34276705000000,
F = 6096922.76188281000000;

function x1(x, y) { return A*x + B*y + C; }

function y1(x, y) { return D*x + E*y + F; }

What more do I have to do to get longitude and latitude? For example C, F represents the top-left corner of the image, right? Then how do I translate C, F to get longitude and latitude?



I want do do it using javascript only, have looked at http://proj4js.org/ but don't really understand what to do with the values I get...




Thursday 27 October 2016

pyqgis - QGIS start-up error message: "Couldn't load plugin MetaSearch due to an error when calling its classFactory() method"


I'm an Ubuntu 16.04 (64-bit) user and I just recently upgraded my QGIS instalation (from 2.14, "Essen", to 2.16, "Nødebo").


I'm not a programmer and I know little about the details of how python works, much less about how QGIS and python relate, so when I get this message at start-up, I don't really know what to do about it:


Couldn't load plugin MetaSearch due to an error when calling its classFactory() method 


ImportError: No module named requests
Traceback (most recent call last):
File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 333, in startPlugin
plugins[packageName] = package.classFactory(iface)
File "/home/jmb/.qgis2/python/plugins/MetaSearch/__init__.py", line 29, in classFactory
from MetaSearch.plugin import MetaSearchPlugin
File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import
mod = _builtin_import(name, globals, locals, fromlist, level)
File "/home/jmb/.qgis2/python/plugins/MetaSearch/plugin.py", line 31, in

from MetaSearch.dialogs.maindialog import MetaSearchDialog
File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import
mod = _builtin_import(name, globals, locals, fromlist, level)
File "/home/jmb/.qgis2/python/plugins/MetaSearch/dialogs/maindialog.py", line 44, in
from owslib.csw import CatalogueServiceWeb
File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import
mod = _builtin_import(name, globals, locals, fromlist, level)
File "/usr/share/qgis/python/owslib/csw.py", line 27, in
from owslib.util import OrderedDict
File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import

mod = _builtin_import(name, globals, locals, fromlist, level)
File "/usr/share/qgis/python/owslib/util.py", line 35, in
import requests
File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import
mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named requests


Python version: 2.7.11+ (default, Apr 17 2016, 14:00:29) [GCC 5.3.1 20160413]
QGIS version: 2.16.1-Nødebo Nødebo, 8545b3b


Python Path:
/usr/share/qgis/python
/home/jmb/.qgis2/python
/home/jmb/.qgis2/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python2.7
/usr/lib/python2.7/plat-x86_64-linux-gnu
/usr/lib/python2.7/lib-tk
/usr/lib/python2.7/lib-old

/usr/lib/python2.7/lib-dynload
/usr/local/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages/PILcompat
/usr/lib/python2.7/dist-packages/gtk-2.0
/home/jmb/.qgis2//python


Note that I can ignore this message and work in QGIS without further issues.




Anyways, can anyone tell me how to fix this?




In case it matters, here are the commands used for upgrading QGIS:


sudo apt-get update
sudo apt install qgis python-qgis qgis-plugin-grass

And these are the lines in my relevant lines from the /etc/apt/sources.list file:


deb http://qgis.org/debian xenial main
deb-src http://qgis.org/debian xenial main





Note2: can this be originated by the non ASCII character in the QGIS version "Nødebo"?




Answer



Had the same error, found the solution here: QGIS issues


Just run:


sudo pip install requests

to install the required module.



GDAL python cut geotiff image with geojson file


I wrote the function which cutting geotiff image by bounding box.


First image is original. At second you can see result off my code. I don't use gdalwarp or any console utilities. But I have no idea how to cut by geojson file. Also I can use only GDAL and numpy modules.


Original image Cutted by bounding box


Cutted by GeoJSON


Here is my code:


import os, sys

from osgeo import gdal,gdalconst,osr


def cut_by_bounding_box(min_x, max_x, min_y, max_y):
xValues = [min_x, max_x]
yValues = [min_y, max_y]
#
# Register Imagine driver and open file
#
driver = gdal.GetDriverByName('GTiff')

driver.Register()
dataset = gdal.Open(filename)
if dataset is None:
print 'Could not open ' + filename
sys.exit(1)
#
# Getting image dimensions
#
cols = dataset.RasterXSize
rows = dataset.RasterYSize

bands = dataset.RasterCount
#
# Getting georeference info
#
transform = dataset.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = -transform[5]
#

# Computing Point1(i1,j1), Point2(i2,j2)
#
i1 = int((xValues[0] - xOrigin) / pixelWidth)
j1 = int((yOrigin - yValues[0]) / pixelHeight)
i2 = int((xValues[1] - xOrigin) / pixelWidth)
j2 = int((yOrigin - yValues[1]) / pixelHeight)

new_cols = i2 - i1 + 1
new_rows = j2 - j1 + 1
#

# Create list to store band data in
#
band_list = []
#
# Read in bands and store all the data in bandList
#
for i in range(bands):
band = dataset.GetRasterBand(i+1) # 1-based index
data = band.ReadAsArray(i1, j1, new_cols, new_rows)
band_list.append(data)


new_x = xOrigin + i1 * pixelWidth
new_y = yOrigin - j1 * pixelHeight
new_transform = (new_x, transform[1], transform[2], new_y, transform[4], transform[5])
#
# Create gtif file
#
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file, new_cols, new_rows, 3, gdal.GDT_Byte)
#

# Writting output raster
#
for j in range(bands):
data = band_list[j]
dst_ds.GetRasterBand(j+1).WriteArray(data)
#
# Setting extension of output raster
#
dst_ds.SetGeoTransform(new_transform)
wkt = dataset.GetProjection()

#
# Setting spatial reference of output raster
#
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#
# Close output raster dataset
#
dataset = None

dst_ds = None


if __name__ == '__main__':
# Imput/output file name and set directory
os.chdir('/home/sant/test/satellite_images')
filename = '20160501.tif'
output_file = '/home/sant/test/20160501_cutted_by_bounding_box.tif'
cut_by_bounding_box(531961.73, 535987.34, 4894164.57, 4888631.61)
print 'cutter.py script done!'


Answer



Here is my own solution. It works for any number of bands, any types of geometry(e.g. multipolygon) and works with images any zones!


import geojson as gj
from osgeo import ogr, osr, gdal

# Enable GDAL/OGR exceptions
gdal.UseExceptions()


# GDAL & OGR memory drivers

GDAL_MEMORY_DRIVER = gdal.GetDriverByName('MEM')
OGR_MEMORY_DRIVER = ogr.GetDriverByName('Memory')


def cut_by_geojson(input_file, output_file, shape_geojson):

# Get coords for bounding box
x, y = zip(*gj.utils.coords(gj.loads(shape_geojson)))
min_x, max_x, min_y, max_y = min(x), max(x), min(y), max(y)


# Open original data as read only
dataset = gdal.Open(input_file, gdal.GA_ReadOnly)

bands = dataset.RasterCount

# Getting georeference info
transform = dataset.GetGeoTransform()
projection = dataset.GetProjection()
xOrigin = transform[0]
yOrigin = transform[3]

pixelWidth = transform[1]
pixelHeight = -transform[5]

# Getting spatial reference of input raster
srs = osr.SpatialReference()
srs.ImportFromWkt(projection)

# WGS84 projection reference
OSR_WGS84_REF = osr.SpatialReference()
OSR_WGS84_REF.ImportFromEPSG(4326)


# OSR transformation
wgs84_to_image_trasformation = osr.CoordinateTransformation(OSR_WGS84_REF,
srs)
XYmin = wgs84_to_image_trasformation.TransformPoint(min_x, max_y)
XYmax = wgs84_to_image_trasformation.TransformPoint(max_x, min_y)

# Computing Point1(i1,j1), Point2(i2,j2)
i1 = int((XYmin[0] - xOrigin) / pixelWidth)
j1 = int((yOrigin - XYmin[1]) / pixelHeight)

i2 = int((XYmax[0] - xOrigin) / pixelWidth)
j2 = int((yOrigin - XYmax[1]) / pixelHeight)
new_cols = i2 - i1 + 1
new_rows = j2 - j1 + 1

# New upper-left X,Y values
new_x = xOrigin + i1 * pixelWidth
new_y = yOrigin - j1 * pixelHeight
new_transform = (new_x, transform[1], transform[2], new_y, transform[4],
transform[5])


wkt_geom = ogr.CreateGeometryFromJson(str(shape_geojson))
wkt_geom.Transform(wgs84_to_image_trasformation)

target_ds = GDAL_MEMORY_DRIVER.Create('', new_cols, new_rows, 1,
gdal.GDT_Byte)
target_ds.SetGeoTransform(new_transform)
target_ds.SetProjection(projection)

# Create a memory layer to rasterize from.

ogr_dataset = OGR_MEMORY_DRIVER.CreateDataSource('shapemask')
ogr_layer = ogr_dataset.CreateLayer('shapemask', srs=srs)
ogr_feature = ogr.Feature(ogr_layer.GetLayerDefn())
ogr_feature.SetGeometryDirectly(ogr.Geometry(wkt=wkt_geom.ExportToWkt()))
ogr_layer.CreateFeature(ogr_feature)

gdal.RasterizeLayer(target_ds, [1], ogr_layer, burn_values=[1],
options=["ALL_TOUCHED=TRUE"])

# Create output file

driver = gdal.GetDriverByName('GTiff')
outds = driver.Create(output_file, new_cols, new_rows, bands,
gdal.GDT_Float32)

# Read in bands and store all the data in bandList
mask_array = target_ds.GetRasterBand(1).ReadAsArray()
band_list = []

for i in range(bands):
band_list.append(dataset.GetRasterBand(i + 1).ReadAsArray(i1, j1,

new_cols, new_rows))

for j in range(bands):
data = np.where(mask_array == 1, band_list[j], mask_array)
outds.GetRasterBand(j + 1).SetNoDataValue(0)
outds.GetRasterBand(j + 1).WriteArray(data)

outds.SetProjection(projection)
outds.SetGeoTransform(new_transform)


target_ds = None
dataset = None
outds = None
ogr_dataset = None

pyqgis 3 - Add Help Menu entry in QGIS 3 from `startup.py`



I would like to add a menu entry in the Help menu pointing to some web ressource, say https://gis.stackexchange.com. The following code executed from the python console works perfect:


from qgis.utils import iface
import webbrowser

def open_gis_se():
webbrowser.open('https://gis.stackexchange.com')

iface.helpMenu().addSeparator()

gis_se_action = QAction('Go to gis.stackexchange')

iface.helpMenu().addAction(gis_se_action)
gis_se_action.triggered.connect(open_gis_se)


Result when typed in the python console:


enter image description here


... but putting it into my startup.py has no effect (Help menu remains 'as it is').


In QGIS 2, the above code put in the startup.py adds the desired menu entry as expected.


Why?



Answer




Great idea


you need place startup.py in C:\Users\\AppData\Roaming\QGIS\QGIS3 and add missing import ,and voilá


from qgis.utils import iface
from PyQt5.QtWidgets import QAction
import webbrowser

def open_gis_se():
webbrowser.open('https://gis.stackexchange.com')

iface.helpMenu().addSeparator()


gis_se_action = QAction('Go to gis.stackexchange')
iface.helpMenu().addAction(gis_se_action)
gis_se_action.triggered.connect(open_gis_se)

enter image description here


Calculating number of pixels in polygon using Google Earthy Engine?


How to calculate the number of pixels in a polygon, especially, when an image contains null-value pixels? I tried the function ee.Reducer.countEvery(), but it does not work.


var geometry = /* color: #d63000 */ee.Geometry.Polygon(
[[[-122.420654296875, 37.96148894979038],
[-122.530517578125, 37.735934739135104],
[-122.244873046875, 37.82276800453933],

[-121.904296875, 37.87481898883636]]]);

// Load a Landsat 8 image.
var image = ee.Image('LANDSAT/LC08/C01/T1/LC08_044034_20140318');


var c=image.reduceRegion({
reducer: ee.Reducer.countEvery(),
geometry: geometry
});


print(c);


Wednesday 26 October 2016

python - Programmitcally Change Text Font for all Labels on MXD


I'd like to know if there is a way to change the text font for all layers and all label classes within a layer on a MXD.


I use a font not often found on all computers for created maps. I love the font, but it becomes an issue when publishing or sharing MXD's. Usually before I share or publish a MXD I need to go through all the layers and set the font to be something found everywhere.


Personally I'd like if the script auto found layers and label classes withing every layer in the TOC. I'm pretty good with python, but I haven't ever found a way to interact with the label's for a layer within ArcPy.


Does anyone have an idea of how to due this?




Answer



You cannot change font style using arcpy, your options are:




  1. using label expressions


    "" & [LABELFIELD] & ""

    "" & [LABELFIELD] & ""



  2. ArcObjects - ITextSymbol Interface




I posted a question on this site before, ArcObjects add halo to label, dealing with labels which will give yoiu a general guild for accessing the ITextSymbol interface. You will just need to change the beginning portion of the code to loop through each layer within the map document and add the pFont.Font property.


arcgis 10.3 - Different watersheds for same area because of "burned streams"


Im trying to analyse where water is flowing which rains down to my project area. This is the DEM I have (my project area is bordered blue, the red marks are my pour points):



enter image description here


I learned from another thread (Could bridges prevent ArcMap from creating correct watersheds?) that culverts could disturb the creating of watersheds. So I used the "DEM Reconditioning" tool of ArcHydro to burn the streams:


enter image description here


Then I used this DEM (called it "agreedem") and the following workflow to generate my watersheds (snap distance set to 3 m):


enter image description here


This was my first result:


enter image description here


Since it is my main aim to analyze what I have to do to keep the water inside the center of my area I started building small dams using the "Build Walls" tool of ArcHydro. After that I rerun my watershed workflow (several times with different dams, this was the best result I could get):


enter image description here


I was pretty happy with the result. Then I asked myself: "Why should I use the "burn the streams" thing in this case?". There are no streams with culverts in this area (though I have the situation with the culverts in other areas I'm looking at) so there should be no reason to do so. I deceided to rerun the watershed workflow for my original DEM (called "dtm_orig"):



enter image description here


As you can see, the watersheds are different from the ones I generated with the "agreedem" DEM. So I asked myself how there could be such a great difference in the watersheds between the ones created basing on the "agreedem" and the ones based on "dtm_orig". I couldn't imagine that bigger and deeper ditches / streams (whats the only difference between "agreedem" and "dtm_orig") lead to this results and analyzed all the steps of my "create watershed" workflow.


I then realized that there are huge differences in the DEMs which are generated in the first step by using the "Fill" tool between the "agreedem"...


enter image description here


...and the "dtm_orig" DEM:


enter image description here


This are the values of one of the points where one could see the difference:


enter image description here


Why is there such a big difference between the both filled DEMs?


As suggested by Martin I ran the Sink Evaluation (on the dtm_orig DEM) today. To be honest I don't know how to interpret the results.



This is the "SinkPoly" Layer:


enter image description here


Same with added "SinkDA" Layer:


enter image description here




coordinate system - Projecting EPSG:4326 data in 2D map?


I like to think I'm relatively well-versed in datums, projections, and coordinate reference systems.



OpenStreetMap data is stored in WGS84 Lat/Lon (EPSG:4326). This CRS is geographic and therefore intended to store locations on a 3D globe. When I see this map what am I looking at?


Looking at Reykjavik it appears that the longitude lines are tending towards the North Pole at the Greenwhich Meridian (or maybe they planned the city in perspective view??) so I assume 0,90 is top-centre of the map.


Does this class as a projection, and if so which one?


Is there one projection that is commonly used for lat/lon data?


I have data stored in Oracle using EPSG:4326 and when rendered by GeoServer without a specified projection it also displays these characteristics.



Answer



OpenLayers uses the term 'EPSG:4326' to mean the Plate Caree projection. Referring to 'WGS84' and EPSG:4326 as a projection has been common for so long that it is a source of confusion. This short-hand has been going on since before Google and OpenLayers came on to the scene. For instance, ESRI have been fudging the terms for as long as I can remember. OpenLayers does not reproject to EPSG:900913 on the fly unless you tell it to, which you have to if you want to mix data with Google base maps because the Google API only uses 900913 (which was invented by Google - the numbers being vaguely reminiscent of the word 'google').


arcgis desktop - Dealing with raster: creating it in projected CRS and exporting it in geographic CRS


When creating a raster from a polygon (the CRS of the polygon is a geographic coordinate system), the raster is created in the projected coordinate reference system of the data frame, which is a projected one. I need to do so because I need to specify the Cell size of the raster in meters.


I also need to deliver this raster (previously converted to polygon) in WGS84 (GCS); I am considering to do it so: Changing the Data frame coordinate reference system to WGS84, and, on the result raster, do Data > Export Data and, for the option “Use the same coordinate system as:”, choose “the data frame”.


Is this an accurate way of dealing with that? I mean, aren’t there modifications in areas or coordinates of the delivered shape layer?





QGIS Server Capabilities Response doesn't contain any Layers



This question is similar as this one QGIS Server not working (problem with fast CGI), but with more details:


I installed the latest QGIS from trunk on Linux (Ubuntu 11.04) in the non-default location '/usr/local/qgis/'. I installed apache2 from the repository.


Next,I followed the steps outlined in this post on Linfiniti's blog (http://linfiniti.com/2010/08/qgis-mapserver-a-wms-server-for-the-masses/), but probably I did not understand the part about telling apache about the library path for my CGI; where do I add the ScriptAlias code to?


I created a qgis project with two vector (shapefiles). Under WMS server tab of the project properties I checked 'service capabilities' and gave the title 'trial'. I copied the qgs project file to the /usr/lib/cgi-bin/VECEA directory.


Trying http://localhost/cgi-bin/VECEA/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities in the browser gets me the 'capability' xml file. It contains the user information, but there is no information on the data layers (the qgis project has two shapefile vector layers) and the bounding box is not defined (all set at 0), is that right (I did define an 'advertised extent' in the qgis project under the WMS tab)


In QGIS I can connect to the server. However, it just gives me one row with the name of the project (ID 0), whereas I would expect it to show the two vector layers that were part of the project. When adding the layer, it is empty.


Edit: and now it stopped working altogether. When trying the getCapabilities, I am getting a 500 Internal Server error. Looking in the apache2 error log I am finding the following:


[Tue Sep 27 20:53:03 2011] [error] [client 127.0.0.1] Premature end of script headers: qgis_mapserv.fcgi /usr/lib/cgi-bin/VECEA/qgis_mapserv.fcgi: error while loading shared libraries: libqgis_core.so.1.8.0: cannot open shared object file: No such file or directory


Edit: Got the 'getCapabilities' part to work again. In http://lists.osgeo.org/pipermail/qgis-user/2011-March/011549.html, it was suggested to add LD_LIBRARY_PATH to /etc/apache2/mods-> enabled/fcgi.conf, adding the line "DefaultInitEnv LD_LIBRARY_PATH /usr/local/qgis_git/lib". Only, I had to edit '/etc/apache2/mods-available/fcgid.conf'


Still stuck with no layers showing when connecting to the server in QGIS





gps - What is the practical limit for lat/long measurement accuracy?


There are several systems for finding lat/long coordinates. This includes GPS satellites, DoD TRANSIT, Russian GLONASS, Galileo, differential GPS, LORAN, land-based systems, and more. Maybe even using laser interferometry directly from Greenwich.


I want to know, using all of this technology, what is the actual maximum accuracy which a location on Earth is measurable?



This is not a question about technology or measurement. The theoretical limit is nanometers because of lasers. But the actual limit will be worse because of the physical characteristics of Earth, it is constantly moving.


So what is is actual best measurement we can make?



Answer



It is currently possible to achieve sub-centimeter accuracy when surveying a point.


To express a latitude and a longitude, it needs to be reference to something. We need a datum: a reference frame and a reference ellipsoid. The reference ellipsoid is simply the mathematical shape that we use to perform the conversion to a geographic coordinate system (lat/long/height). The reference frame is the set of physical points (stations) that serve to locate and orient a coordinate system. For example, the ITRF2014 (International Terrestrial Reference Frame) is defined with the coordinates of hundreds of stations around the world, as well as their velocities (because of crustal movements). The ITRF2014 is globally non-rotating relative to the crust, is centered on the geocenter, and has no translation. It is an ECEF (Earth-centered, Earth-fixed) frame of reference.


This most recent realization of the ITRS (ITRF2014), with modern geodetic techniques (mainly VLBI, DORIS, SLR and GNSS) has a very good overall accuracy and consistency, and the errors are typically only a few millimeters. Countries often use other local or national datums to express positions, using reference frames that locally follow the crust and reference ellipsoids that better approximate the shape of the Earth in that country. Accuracy and consistency of these datums can vary a bit depending on which techniques were used to realize their reference frames.


Now for the practical measurement part, if you wish to determine the lat/long of a point, RTK GPS, and even traditional surveying techniques can be used to achieve a reasonable accuracy of a few centimeters pretty quickly. These techniques require a nearby known reference point or broadcasting station for GPS, usually from a densification network, a local geodetic network or some other previously measured point. With more advanced GPS survey, if you remain static for several hours or days and collect enough data from satellites and post-process your data using nearby ground station corrections, it is generally possible to get down to a few millimeters in accuracy. Sub-centimeter accuracy surveys can help analyzing crustal deformation, tectonic motion, isostatic rebound, sea level rise, etc.


arcgis desktop - How do I combine bands in a tiff then convert to grayscale?


I have a large mosaic of tiffs in which I need to do the following: -Combine bands 4,3,2, for an infrared image then convert the image to grayscale without producing another dataset -This is easy on QGIS but I need to execute this on ArcGIS 10.1


Is there a way to render the image like on QGIS?



Answer




There is no direct method, but there are two ways of achieving this.


First, it can be done using the image analysis window (Windows > Image analysis). Select your image, then apply your function (fct button in "Processing") on it. To select a function (in your case, "composite band" and/or " grayscale"), right-click on the identity function and select "insert"). Mainly interesting if you have many functions to apply.


enter image description here


On the image, you see my "on the fly" grayscale image from a multispectral image.


The other method consists in computing your grayscale image (band1 + band 2+ band3)/3 in raster calculator with an in memory output. No need to create a composite raster before.



in_memory/grayscale



As a remark, slightly out of ArcGIS, ArcGIS supports the vrt format. So you could also create a vrt and open it in ArcGIS.


python - Specifying region of interest for google earth engine to calculate Jeffries-Matusita separability?


My objective is to compute Jeffries-Matusita separability using google earth engine python api. I have never worked with ee before, so I am trying to follow this github.


In it, to import roi it says:


table = ee.FeatureCollection('users/mortcanty/supervisedclassification/train')

trainData = image.sampleRegions(table,['CLASS_ID'])

What surprised me here is that the file here has no extension. So I tried to use shp (and only shp,not its supporting files).


However, subsequently when I try to parse 'table' in this function


def jmsep(class1,class2,image,table):
# Jeffries-Matusita separability
table1 = table.filter(
ee.Filter.eq('CLASS_ID',str(class1-1)))
m1 = image.reduceRegion(ee.Reducer.mean(),table1)\
.toArray()

s1 = image.toArray() \
.reduceRegion(ee.Reducer.covariance(),table1)\
.toArray()
table2 = table.filter(
ee.Filter.eq('CLASS_ID',str(class2-1)))
m2 = image.reduceRegion(ee.Reducer.mean(),table2)\
.toArray()
s2 = image.toArray() \
.reduceRegion(ee.Reducer.covariance(),table2,15)\
.toArray()

m12 = m1.subtract(m2)
m12 = ee.Array([m12.toList()]) # makes 2D matrix
s12i = s1.add(s2).divide(2).matrixInverse()
# first term in Bhattacharyya distance
B1 = m12.matrixMultiply(
s12i.matrixMultiply(m12.matrixTranspose())) \
.divide(8)
ds1 = s1.matrixDeterminant()
ds2 = s2.matrixDeterminant()
ds12 = s1.add(s2).matrixDeterminant()

# second term
B2 = ds12.divide(2).divide(ds1.multiply(ds2).sqrt())\
.log().divide(2)
B = ee.Number(B1.add(B2).project([0]).toList().get(0))
# J-M separability
return ee.Number(1).subtract(ee.Number(1) \
.divide(B.exp())) \
.multiply(2)



a = jmsep(5,9,image,table).getInfo()

Gives error:



Collection asset 'C:\Users\train.shp' not found



I suspect this is due to 'shp' file not being appropriate.



Answer



The problem is that the table needs to be on the server side for Earth Engine to be able to use it.


The step prior to this script is to add the .shp as an asset in your Earth Engine user.



If you look at the "Assets" tab, there is a red NEW button. Here you are able to load your shapefile as a table asset that you can use later in your script.


shapefile - Given a collection of polygons, how do you find which one contains a given point? Java or Python


I'm looking for Java or Python libraries that will help me answer this question efficiently: given a collection of polygons (via shape files or GeoJSON), which one of them contains a given point? The overall objective is to compute counts of points by polygons.



Answer



Use can do this in Java with this library: https://github.com/Esri/geometry-api-java


OperatorDisjoint with polygon vs point let's you answer the question for a point vs polygon.


To make it more efficient you need to make a spatial index that would allow to limit the number of points vs the number of polygons before testing with OperatorDisjoint.



If you have enough memory, you could use QuadTree class as the spatial index, storing either all points, or all polygons (performance depends on data configuration).


Say you put points in the QuadTree, then for each polygon, query the QuadTree for points that intersect the extent of the polygon, then use OperatorDisjoint on each point/polygon. If there are more than 10 points to test against given polygon, accelerate the polygon, using OperatorDisjoint.accelerateGeometry, before running the OpeatorDisjoint.execute. Accelerating builds an index on the polygon itself, making point in polygon test much faster. If polygons stay in memory, call OperatorDisjoint.deaccelerateGeometry on the polygon, once done with it.


If you query points and polygons from a database, then it probably already has spatial index, and you can skip QuadTree part.


The library also provides import from GeoJSON and Esri shape.


arcgis 10.1 - Using ArcPy to Switch Layout Templates?



Using ArcMap 10.1, since there is not a way to control portrait or landscape of the layout view through python, is there a way to switch layout templates without opening a new mxd through python?


Since there is a change layout template button it would be nice to have a python function to do the same.




arcobjects - Programmatically edit/update metadata in ArcGIS?


Has anyone succeeded in programmatically updating metadata in ArcGIS 10? Considering using Python/arcpy but ArcObjects (C# or Python/comtypes) is also a possibility.



I need to update both the FGDC and the ArcGIS-ISO format metadata, and whatever solution is used needs to be able to retain the existing (non-blank) elements along with the added elements, except where they are in conflict in which case the added elements overwrite the existing elements.



Answer



The easiest way to do this from arcpy is to create an xml file using python and then invoking Import Metadata (Conversion). However, this will overwrite everything.


An alternative is to use ArcObjects to obtain an IName to the dataset, cast to IMetadata, and edit the IPropertyset.


arcgis 10.0 - Longitudinal Issues with NetCDF data


I am seeking advice on working with world data sets where the GCS has a positive range (0 degrees to 360 degrees). The data I am working with are NetCDF oceanographic data and have positive coordinate values as mentioned. When displayed in a normal GCS WGS84 in ESRI's ArcGIS 10, it is offset from other data that exist in the range of -180 to 180. If I reproject it, gaps occur around the prime meridian (mostly due to its origins, sometimes as far off as 25 degrees W). My thought is to create a custom projection/coordinate system that has the positive values and hope that other world data sets reproject on the fly just fine. Does anyone see any issues with this or have other solutions? Hope this makes sense. I appreciate any help in advance.


** Update ** This is a screenshot of a countries data set with graticule on top to show the normal CRS of GCS. Just importing the NetCDF data as-is brings it in but the values don't start until 20.5 degrees east of the Prime Meridian.



Normal GCS http://grafa.co/rnd/img/GCS_Normal.png


If I choose to use the NetCDF's CRS (which is actually the same) it will reproject the world data on the fly if it is defined. Note the graticule is undefined so it does not reproject. Sort of reprojected GCS http://grafa.co/rnd/img/GCS_Other.png


But then if I try to reproject everything to a normal GCS with negative values, it's like wrapping a flat map of the data around the world and then it disappears when it hits the Prime Meridian. Sort of reprojected GCS http://grafa.co/rnd/img/reprojected.png


Now, I know there are no values from 0-20.5 as is stated in the metadata. But why can't the values from rest of the data display in the gap? I even tried a Shift in the Raster Tools to no avail.




arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...