texts: Calculating ratio between "rows" in ArcGIS field calculator

texts:
I have a problem similar to Dividing two or more rows in single attribute table? but that is for QGIS whereas I have to perform it in ArcGIS.
I have to calculate the Ratio between rows (sequentially, i.e., row1/row2, row2/row3 and so on) of a specific field based on the ID field and have to assign the results to a separate field.
Here is the code and the error message:
 def CalcRatio(ID,Frequency):
    vDict = {}
    with arcpy.da.SearchCursor("Test",['Frequency'],'ID = {}'.format(ID)) as sCur:
    for Row in sCur:
    vDict[Row[0]]=Row[1] # Dict[id] =Frequency
    try: value = ( vDict[1] / vDict[2])
    except: value = -1
    return value

Message:
ERROR 000539: Error running expression: CalcRatio( 1 , 50680) Traceback (most recent call last):
File "<expression>", line 1,
in <module> File "<string>", line 5,
in CalcRatio IndexError: tuple index out of range

Here is the data
OID Id  Frequency
0   1   50680
1   2   49740
2   3   48620
3   4   48300
4   5   48004
5   6   45600
6   7   44384
7   8   42720
8   9   41890
9   10  39700
10  11  38530
11  12  38207
12  13  38106
13  14  37303
14  15  37205
15  16  36527
16  17  35301
17  18  31674
18  19  28890
19  20  24612


texts: Create TIN from Elevation Points in Python

texts:
I want create a TIN from point features with elevation using a python script.

This is my script:

import arcpy
from arcpy import env
import json as simplejson
import urllib
import sys,os

arcpy.CheckOutExtension("3D")

env.workspace = r"D:\test\Pythontest"

OutputTin = "PointFTin"
Delaunay = "constrained_delaunay"

arcpy.CreateTin_3d(OutputTin,{},"D:\test\Pythontest\PlotGrid_label_prj.shp"+ "Shape.Ele masspoints", Delaunay)

When I run this script I get this error:

RuntimeError: Object: Error in executing tool

Where is it going wrong?


texts: In GeoWebCache is there a way to disable the geowebcache.xml backups? (.bak)

texts:
I have a very large geowebcache.xml file that is edited constantly. Every time the system makes a change to geowebcache.xml via the REST API GWC creates a backup of the geowebcache.xml file. These backups are taking up massive amounts of space each day and aren't being rotated by GWC.

I've checked the code and it seems it's only supposed to keep 10 backups, but that isn't happening. Is there any configuration setting to disable these backups?

The alternative is just to create a cron-job to periodically remove them I suppose. P.S we're using standalone GWC 1.8.


texts: clip low resolution raster with polygon keeping exact boundary

texts:
For visual purposes (not analysis) I want to "clip" a low resolution raster using a shapefile with the exact boundary of the shapefile. ST_Clip doesn't seem to allow cutting across the inside of a pixel so it gives jagged edges since whole pixels are preserved. Is it possible to do this type of clipping without resampling to a higher resolution?

Clip maybe isn't the correct term to use which might be why I'm having trouble finding an answer. Something like a cookie cutter effect. In the end I just want to serve out a png of the bordered area.


texts: Apply function to calculate NDVI on list of images

texts:
I'm looking for a simple for-loop to apply a function to create NDVI images from an original list and save those in a new folder.
I've created a function and show how I would apply it on one image but I would like to loop through a file that contains multiple images AND save those to a new folder.
image_file = "Path/to/image.tif"

raster = rasterio.open(image_file)

red = raster.read(3)
red_edge = raster.read(4)
nir = raster.read(5)

np.seterr(divide='ignore', invalid='ignore')

def calc_ndvi(nir,red):
    '''Calculate NDVI from integer arrays'''
    nir = nir.astype('f4')
    red = red.astype('f4')
    ndvi = (nir - red) / (nir + red)
    return ndvi

ndvi = calc_ndvi(nir,red)


texts: Styling river charts with TileMill - issue with filling a polygon made of disjoint line segments

texts:
I am trying to style inland river charts produced by the Army Corps of Engineers. I imported the vectors into Postgres, and I’m using TileMill to style the data.

I am having trouble correctly filling the inner contour for the river, because the contour is made of line segments that are broken by features like bridges, and the contour fill spills out.

How can I fill the contour properly and avoid this issue? Some sort of processing in Postgres? Is this a solved problem… a good way to make gapped line segments into polygons for filling?

You can see a screenshot of the issue here:


texts: Exporting set of lat/long coordinates to SVG using tool or code?

texts:
Suppose I have a set of coordinates in latitude and longitude.

What would be the easiest way to export those to an SVG file?

I'd rather do this programatically (any language or library will do) rather than manually, but can revert to that if necessary.


texts: How do I convert a set of latitude/longitude coordinates with timestamps to a GPX tracklog?

texts:
Converting tracks from a GPX track to a CSV table or some other format is relatively easy, but I'm looking to go the other way. I have a table of point locations logged with a handheld GPS, with timestamps. I need to convert them to a GPX file with a tracklog containing those points so they can be used with some other software to geotag some photos. I can convert them to GPX waypoints easily enough using QGIS/GDAL, but not to a GPX tracklog, the timestamps are lost.

Here's the header of my CSV:

Lon,Lat,Altitude,NewTrack,TimeStamp,Datum,InternalTime


texts: My QGIS (2.8.2) joins always result in an outer join.  How do I create an inner join with QGIS?


texts: Spotty and Impossible Route (Closest Facility)- ArcGIS Network Analyst

texts:
I am trying to find the closest facility for parcels (converted my polygons to points), to existing sewer lines. In the image you can see that some sections of the green routes are just very short segments, that never reach the existing sewer lines (I converted these to points as well for use in Load Locations). How can this be a valid route when it doesn't connect the parcels to the existing sewer?

I have set up topology on the proposed sewer lines to make sure they are all connected.

I have also tried using the Advanced Options in Load Locations (Incidents which is Parcels), and setting my 'Snap to Position Along Network" with an offset of 20 meters since that is common for parcel points to be 20 meters from proposed sewer lines. This didn't change anything.

I am using Use Feature Geometry option in Location Position, and leaving default value of 5000 meters.

This is related to my previous question:
Calculate distance from polygon to nearest existing line in geometric network, along network routes


texts: How to separate shapefile by attributes into a bunch of new shapefiles?

texts: I am trying to work with census block data at the moment, and I want each of the 4178 blocks making the entire state to be their own shapefiles. How can I do this without having to manually query or iterate through each feature and turn each census block into a shapefile with the "Extract Data" function?


texts: Can you reorder a classified vector layer in QGis?

texts: I have a GML file that I'm displaying in QGis - it contains a set of contours and they are written out in such a way that the lowest contour (largest polygon) is drawn last thus covering up all the smaller polygons. While I can work round this by writing out the GML in reverse order I was wondering if there was a way to reverse the draw order directly inside QGis.


texts: QGIS-I have 37 polygons, can I bulk create 2 other sets of 37 that are 1 px larger all the way round and 1px smaller

texts:
I have 37 polygons covering a Lakes spatial extent at a particular point in time. I need to bulk create 2 sets of polygons. One 1 px larger all the way around and one 1px smaller to calculate the area so I can work out the error for the initial polygons created.

Other suggestions about working out the error are welcome!


texts:
I have a DXF file which is not geocoded, but it would displayed correctly in QGIS but not at the right place.

There are a hugh amount of layers in subgroups.

Selecting one by one would take hours.

How can i move it to the right point / coordinate?

I do not have AutoCAD.

Is it possible to achive it with QGIS?


texts: Merging shapefiles and creating unique feature wherever shapefiles' features overlap using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.5.

I have a shapefile for impervious surfaces and a shapefile for tree canopy. I need to separate the tree canopy features into tree canopy that overlaps (is assumed to be above) impervious surfaces vs. tree canopy that does not overlap any impervious surface, and specifically I need to be able to calculate the area of each. Ideally I'd also merge these layers together to get a single shapefile that shows tree canopy over pervious cover, tree canopy over impervious cover, and impervious cover with no tree canopy above it, with the ability to calculate area for each class.

Here's a screenshot of what I have so far:

How can these shapefiles be merged so to create a unique feature or attribute wherever the original shapefiles' features overlap?


texts: Generating multiple water level scenarios with only DEM data in QGIS

texts:
I want to depict certain water level scenarios (from -0,40mNN to +0,20mNN) by using a two dimensional hypsometric map. Only provided data are DEM-files of the AOI.
I managed to create the hypsometric map, but how do I do the water level scenarios with only the DEM data?
I am using the desktop version of QGis (Version 3.10.5). Here is a picture of what the map should look like as a result:


texts: turf.intersect of a turf.polygon and a turf.linestring?

texts:
With the JavaScript library turf, when I try to use turf.intersect(polygon, linestring)) I get the following error:

ReferenceError: 'coords' is undefined
at jsts.geomgraph.GeometryGraph.prototype.addLineString (http://localhost:6235/Scripts/turf.min.js:15:12334)
at jsts.geomgraph.GeometryGraph.prototype.add (http://localhost:6235/Scripts/turf.min.js:15:11273)
at jsts.geomgraph.GeometryGraph (http://localhost:6235/Scripts/turf.min.js:15:9430)
at jsts.operation.GeometryGraphOperation (http://localhost:6235/Scripts/turf.min.js:7:9064)
at jsts.operation.overlay.OverlayOp (http://localhost:6235/Scripts/turf.min.js:7:16807)
at jsts.operation.overlay.OverlayOp.overlayOp (http://localhost:6235/Scripts/turf.min.js:7:17248)
at n.prototype.getResultGeometry (http://localhost:6235/Scripts/turf.min.js:13:17602)
at n.overlayOp (http://localhost:6235/Scripts/turf.min.js:13:17228)
at jsts.geom.Geometry.prototype.intersection (http://localhost:6235/Scripts/turf.min.js:2:20228)
at e.exports (http://localhost:6235/Scripts/turf.min.js:17:18325)

What am I missing?


texts: Why is {maxPixels : } set to 1024 for MODIS and 80 for Landsat?

texts:
I am currently trying to understand how images are resampled and reprojected in google earth engine and I came across this tutorial: https://developers.google.com/earth-engine/resample

Why is maxPixels:1024 and maxPixels:80 for MODIS and Landsat respectively? I searched specifications but did not find this information.


texts:
I have two layers that look like this:

I want the intersection (or whatever the correct term may be) to look like this where the output polygons are the same as the original base layer, just with added attributes if there is any kind of intersection. Preferably, there'd be a way to say something like, if there is over 50% coverage, take the property, if not, no property.

Similarly, this area has a lot of intersections where it shouldn't. (The C4 applies to the dark shapes and the numbers apply to the lighter shapes. The darker wide line is where the two number polygons overlap.)

Ideally, the "fuzzy" intersection would work similarly here where the polygon below gets the added property of the polygon above that overlaps the most with it. Again, the output should be with the same polygon shapes as the base layer.

How can I achieve these results? I've looked into Union, Intersect, Clip, Identity and such, but none seem to be right. Perhaps there is an option in one of these that I am missing? I can use either ArcGIS or QGIS.


texts: How do I export training samples from IDRISI?

texts:
I have been tying to implement this paper on an SVM classifier wherein I need to extract training samples from an image and run an SVM classifier on it. The paper says the following:

  IDRISI Andes was used for the preliminary processing such as
  identification of  training samples. These data were then exported to
  MATLAB (version 7) for SVM classification.

However, I cannot figure out how to export the training sample data in a format readable by MATLAB (or even python, for that matter). The help page I found doesn't seem very useful.

I found answers here that tell me eCognition has a utility to do this, and so does grass GIS. Is it recommended that I use one of these softwares, or is there any way to do it in IDRISI?

(I also found this question which does not have any answers.)


texts: How to convert from EPSG:4326 to WGS84/UTM 32U?

texts:
I am unable to convert from the project's coordinate system to a different one using QGis 1.8.0.

I have a .shp file with a certain coordinate system information (WGS84 EPSG:4326 according to qgis). I wish to convert this to WGS84/UTM 32U.
Now, WGS84/UTM 32N but not [...] 32U exists in the KBS list if I use [...]32 N it does not change anything.

If I compare the coordinates in my project in WGS84 with that from 2 different webmap applications in WGS84/UTM 32U there are these following difference:
3607000 E / 5738000 N in the .shp and
607000 E / 5738000 N in the .shp in the web application with an seemingly difference of ~100m west and a ~2000m north shift of the .shp.

Does someone have an idea what the problem might be?
Can I download more Coordinatesystems (i.e. WGS84/UTM 32U) for QGIS somewhere?

Thank you very much for any help.


texts: How to use ground overlay google map android for polygon area to fit the map image in it?

texts:
I have a park area that is not rectangle or a simple shape. It is just like polygon area. I could load sample image on it for testing, in my android app.
There it is :

mGroundOverlayRotated = googleMap.addGroundOverlay(
    GroundOverlayOptions()
        .image(mImages[0])
        .position(LatLng_NEAR_NEWARK, 255f, 84f)
        .bearing(110f)

It is worked for simple shape like rectangle.

LatLng_NEAR_NEWARK is latitude/longitude of center of park.
255f is width of area in meter.
84f is height of area.

I cannot using width/height for the area like this.
I can't figure out how to use Position function in my code to fit it.

I need something like this :

https://stackoverflow.com/questions/26250481/show-only-ground-overlays-map-android

I saw something about KML format. I'm not sure about it. Should I use SVG/PNG format for this area instead of JPG?

What approach do I use?


texts: Import CSV and join with another table to add geometry in QGIS

texts:
This relatively simple task has me stumped in QGIS so want to check if I'm missing something.

We get a lot of data supplied in spreadsheets. These datasets have no spatial information but will contain a field we can join with a spatial dataset (e.g postcode, hospital ID, census code etc)

In Mapinfo we have a few options:

Open CSV (or XLSX), open the spatial TAB and use the "Geocode" function to spatialise the CSV based on the relevant ID column.
Use the "SQL Select" function to join the CSV to TAB and then save the query results as a new TAB.
If the CSV has coordinates (rarely) or we have joined with a TAB that does contain coordinates, but no geometry (more likely) then we just use "Create Points", select the XY columns and we're done.

Either of those options takes literally seconds to do and this is something we do a lot.

My current workflow(s) in QGIS is like this, which seems long winded/doesn't produce satisfactory results:

Open CSV or XLSX with no geometry (e.g List of postcodes with a number we need to theme)
Open postcode layer with ~2.6 million points
Expose the XY coordinates as fields on the postcode layer
Create join on the CSV layer and bring across the XY fields
Save the CSV join as a CSV
Import new CSV this time selecting the new XY fields to create points

Another way I have tried in QGIS is:

Open CSV and postcode table
Create join on postcode table and bring across the field we need to theme, lets call it "Total Visits".
We now have a table of 2.6 million irrelevant postcodes with only a handful of records containing values for "Total Visits" (i.e the ones we want to map)
I try and filter the postcode table down to only show the records with a "Total Visits" value but the joined fields are not available to the filter tool.
Save the full joined postcode layer to a new file
Open the new table and filter to only show relevant postcodes
Save filtered results to new file

I also tried "Join Attributes Table" from the Processing toolbox - while this was a bit quicker than the second option above, you still end up with the full postcode table that needs to be filtered down to only contain the joined records.

Also tried the MMQGIS plugin "Join By Attribute" but when run it claims my postcode layer has no geometry... but it does.

None of those options are very efficient or produce the desired output. Please tell me I'm missing something simple like a hidden option or a plugin (I have searched for various terms but not found a suitable plugin that works in this particular scenario).

At this stage it seems quicker to:

A) Just do it in Mapinfo beforehand
B) Load the CSV into a SQL database and do the join there

So is there a way to essentially do a full outer (or left outer) SQL join within QGIS?


texts: OpenLayers2 addUniqueValueRules default style not displaying

texts:
I'm trying to alter the colour of displayed points using addUniqueValueRules in OpenLayers. If I don't use addUniqueValueRules all my points display the default colour (green). However, when I add addUniqueValueRules, only those points that meet the rules are drawn.

The relevant code:

var styleMap = new OpenLayers.StyleMap({
    fillOpacity: 1,
    fillColor: "green",
    pointRadius: 5
});

var styleLookup = {
    "Hög": { fillColor: "red" },
    "Måttlig": { fillColor: "yellow" },
    "Undersöka": { fillColor: "cyan" }
};

styleMap.addUniqueValueRules( "default", "risk", styleLookup );

wfs = new OpenLayers.Layer.Vector("Editable Features", {
    strategies : [new OpenLayers.Strategy.Fixed(), saveStrategy],
    projection: new OpenLayers.Projection("EPSG:4326"),
    protocol: new OpenLayers.Protocol.WFS({
        version: "1.1.0",
        srsName: "EPSG:4326",
        url: "http://localhost/cgi-bin/tinyows",
        featurePrefix: "tows",
        featureNS :  "http://www.tinyows.org/",
        featureType: "tradcopy",
        geometryName: "geom"
    }),
    styleMap: styleMap
});

Can anyone help me draw all the points - both default and the ones that meet the rules?


texts: Polygonized a land use map (raster) in QGIS but the output is just a black shapefile. How do I fix this?

texts:
So, I have a land use map that I need to digitize. The changes in the map is very fine so it will take a long time to manually digitize it. When I used the Raster > Conversion > Polygonize (Raster to Vector) function, the output shapefile is just black. I tried bands 1 to 4 but same result.

I also tried SAGA's Raster values to points and GRASS' r.to.vect but the result is the same: a black shapefile.

Can anyone help?


texts: Why does ArcMap show quantities for values that are not in attribute table?

texts:
I've joined a csv file of values (0-100) (created in Excel), with local authority boundaries for the UK in ArcMap 10.1. I want to display the percentage of these areas not receiving 3G mobile coverage.

I've gone into properties and then selected quantities. I've selected the data I want to display and it displays huge values which are not in either the original data set or in the attributes table. Why is this? Here is a screen shot (note the range):

When I go into the attributes table for the selected variable ('t') and then sort descending, the largest value is 99. Why does it think there are values up to 1965621253?

I've not managed to find an answer online but so far I've:

Checked that Excel isn't adding on extra values (I copied and pasted just the values numerous times)
Generated a new csv file
Tried to join the data numerous times to see if this made a difference.


texts: Using Python to input datasource of layer into table field

texts:
I am using a bit of Python code to create an attribute table field and then I want to populate that field with the layer's datasource.

I have been using addField management and then calculateField management but I don't know how to reference the datasource so I can have it entered into this field.

This is the key bits of the expression:

inFeatures = "feature" #references the feature in the mxd
inFieldName = "fieldname" #references the field in the feature attribute table
expression = ***this is where I am stuck***
arcpy.CalculateField_management(inFeatures, inFieldName, expression)

It works if I enter a literal string as the expression but I want it to reference the datasource for the feature.


texts: Getting Script Error has occurred in script on this page when running polygon to raster conversion?

texts:
i keep getting the script error problem every time i try using polygon to raster tool in ArcGIS 10.1. Please help me how I can solve this?

screenshot embedded


texts: R - how to create a pre-determined number of identical square polygons to use for a choropleth map?

texts:
I am trying to plot crime data density (again!) over a city map, say of NY.
As a well known problem there are plenty of examples on this (http://www.obscureanalytics.com/2012/12/07/visualizing-baltimore-with-r-and-ggplot2-crime-data/).  These methods plot the crime density through isoclines, while I need to represent it through identical density squares of a pre-determined area (and the area / side length may change from one iteration to the other). The reason for representing crime density through squares is that the square are the hotspot areas to be patrolled. The size will influence the overall amount of police people required.

This is what I am trying to achieve:

I would like to be able to create identical square polygons with pre-determined size to overimpose to the map (is it a raster? apologies but I've just started to learn to spell GIS!)
I would like to use the above squares as items to colour as in a choropleth map (i.e. different colouring in relation to frequency of crime in the area), probably using ggplot2 or similar.
This should allow me to see how the density of crimes per square kilometre varies changing the size (i.e. the area) of the square.

I do not have a clue if it is possible to use R to create regularly shaped squares polygons of a pre-defined size to use for this (as the code snipped below attests). Any help or links to examples are welcome.

I would be glad to get some indication on alternative ways to calculate the density.  I have used the stat_density2 (part of ggplot2) but maybe there are better / faster ways?

This is where I got to:

library(rgdal)
library(raster)
library(sp)
#NY boroughs shapefile downloaded from NY website
shp <- readOGR(dsn = "nybb_14a_av", layer = "nybb")
r <- raster(extent(shp))
res(r)=0.05
# using BoroCode as an experiment...
r <- rasterize(shp, field="BoroCode", r)
plot(r)
plot(shp,lwd=10,add=TRUE)
#don't know the result of the above: the laptop basically hangs processing
#plot(r) :)


texts: Optimal path proposed by pgr_astar is not the optimal one?

texts:
I'm using the function pgr_astar function of PostgreSQL to find the optimal path between two points in the dataset of Chicago.

The query I'm using is:

SELECT a.seq, a.id1 AS node, a.id2 AS edge, b.source, b.target, b.cost, b.reverse_cost, b.km, b.estimated_time as time, b.ang_elev, b.geom_way
 FROM pgr_astar
('SELECT id, source, target, cost as cost, x1, y1, x2, y2, reverse_cost FROM chicago_2po_4pgr', 48,89654, true, true)as a LEFT JOIN chicago_2po_4pgr as b
ON a.id2 = b.id
order by seq

I want to find the most fuel efficient path, the columns cost and reverse cost present the fuel consumed by each segment.
The path proposed by the function is the green one.
BUT, I found a path that have a total fuel consumption less than the proposed one (the yellow path)

Any explanation please? Why didn't the function propose me the yellow path?


texts: How do I draw smooth vertexes around curves in ArcMap 10+?

texts:
I have some polygon layers (parcels) that have vertexes drawn around curves, as if traced.
How do I click and drag to draw vertexes in a smooth shape? Big difference between what was originally there, and what I drew in. Hopefully this helps illustrate what I am wanting to accomplish.


texts: How to add image ol.control.ZoomToExtent in openlayers 3?

texts: In ol.control.ZoomToExtent,Openlayers3 use 'E' character for view (http://openlayers.org/en/v3.5.0/examples/navigation-controls.html). But i want to add image.How to do this ?


texts: Measure shortest distance between raster cell to raster cell of another raster in R

texts:
I would like to measure the shortest distance between a raster cell center and the cell center of another raster.

The low resolution raster (dark grey in the picture) has NA cell values and I would like to assign to it the distance in meter to the closest raster cell of the small raster (the one with the small cells in light grey).

I tried to use the QGIS Proximity(raster distance) tool but somehow it is not working as I wish.

How do I assign to the big raster cells the distance to the nearest raster cell of the small raster?

There are more than one million raster cells which excludes the possibility to vectorize the raster cells (not enough computing power).

Does anyone know an R approach?


texts: Most suitable interpolation method - regridding

texts: I'm looking for guidance on the choice of interpolation method. The data I want to interpolate is air quality data, but not from stations - it is a kind of a model in a 1x1 km grid (where each point represents the average value of the respective cell) and I want to interpolate it to 50x50 m (in the fine grid, each point should represent the value at exactly that location). I tried to do it in R using bilinear (fields::interp.surface) and bicubic (akima::bicubic) methods and I've also considered kriging - but I have a feeling that that would not be suitable for regridding. There are some (not so significant) differences in results, but I would really like to know which approach is the most correct. Maybe even another method? There is just one condition, that I have to be able to do it in R (just as the rest of my code is).


texts: QGIS 2.4 mapRenderer.extent() does not equal mapRenderer.setExtent()

texts:
Could anyone help explain this QGIS 2.4 python console code? I would expect the return value of .extent() to be the same as the bounding box given to .setExtent():

mapRenderer = iface.mapCanvas().mapRenderer()
llx=729200
lly=4361200
urx=736400
ury=4368400
rval=mapRenderer.setExtent(QgsRectangle(llx,lly,urx,ury))
print "  lower left = "+str(llx)+","+str(lly)
print "  upper right = "+str(urx)+","+str(ury)
print "rval:"+str(rval)
print "mapRenderer.extent():"
e= mapRenderer.extent()
print e.toString()

>>>execfile(u'C:/Users/Tom/Documents/extent_test.py'.encode('mbcs'))

  lower left = 729200,4361200
  upper right = 736400,4368400
rval:True
mapRenderer.extent():
727776.0840998685453087,4361200.0000000000000000 : 737823.9159001314546913,4368400.0000000000000000


texts: How to snap to a geojson-vt layer (Vector Tile Layer)?

texts:
I'm using geojson-vt to render a large geojson (originated from OSM) file. This file contains all road from a state (extracted with overpass-turbo way[highway="*"]). I'm buinding a domain-specific application that needs to enable the user to add new roads to the existing geojson. However, I'm unable to snap the drawing cursor to the geojson-vt layer.

To illustrate that, consider a simplified demo with this issue (see fiddle below). This demo includes the geojson of countries limitsand add a drawing cursor that should snap to their boundaries.

JSFiddle: http://jsfiddle.net/mju68Lr4/

  var url = 'https://openlayers.org/en/v3.20.1/examples/data/geojson/countries.geojson';
    fetch(url).then(function(response) {
      return response.json();
    }).then(function(json) {
      var tileIndex = geojsonvt(json, {
        extent: 4096,
        debug: 1
      });
      var vectorSource = new ol.source.VectorTile({
        format: new ol.format.GeoJSON(),
        tileGrid: ol.tilegrid.createXYZ(),
        tilePixelRatio: 16,
        tileLoadFunction: function(tile, tileCoord) {
          var format = tile.getFormat();
          var data = tileIndex.getTile.apply(tileIndex, tileCoord);

          var features = format.readFeatures(
            JSON.stringify({
              type: 'FeatureCollection',
              features: data.features
            }, replacer));
          tile.setLoader(function() {
            tile.setFeatures(features);
            tile.setProjection(tilePixels);
          });
        },
        tileUrlFunction: function(tileCoord) {
          return [tileCoord[0], tileCoord[1], -tileCoord[2] - 1];
        }
      });
      var vectorLayer = new ol.layer.VectorTile({
        source: vectorSource
      });
      map.addLayer(vectorLayer);
      var snap = new ol.interaction.Snap({
          source: vectorSource
      });
      map.addInteraction(snap);

    });

I expected that after adding the geojson-vt (tile layer) I would be able to snap to it, as follow (inside the callback):

  var snap = new ol.interaction.Snap({
      source: vectorSource
  });
  map.addInteraction(snap);

However, I'm getting the following error output (in chrome console):

Uncaught (in promise) TypeError: this.source_.getFeatures is not a function
    at ol.interaction.Snap.getFeatures_ (ol-debug.js:70564)
    at ol.interaction.Snap.setMap (ol-debug.js:70669)
    at ol.Map.<anonymous> (ol-debug.js:41259)
    at ol.Collection.boundListener (ol-debug.js:3441)
    at ol.Collection.ol.events.EventTarget.dispatchEvent (ol-debug.js:3859)
    at ol.Collection.insertAt (ol-debug.js:12466)
    at ol.Collection.push (ol-debug.js:12490)
    at ol.Map.addInteraction (ol-debug.js:41313)
    at (index):106

Am I missing something? It seems that snap is not being able to handle the vector tile layer elements since they are "dynamic". Hence, is it possible to fix this issue and snap the drawing cursor to the geojson-vt layer (a tile layer)? If so, how can I do that?


texts: How do I use variables inside a arcpy toolbox function?

texts:
I have a point file and I need to make a lot of buffers in specific increments but I can't figure out how to interate through a list inside the buffer function. Here is what I have:

import arcpy
city_pts = "C:\Urban_Land_Project\Cities\WI_cities.shp"
buffer_sizes = [1,2,4,7,10,15,20,30]
for x in buffer_sizes:
    arcpy.Buffer_analysis(city_pts, "C:\Urban_Land_Project\City_buffers\WI_city_buffers*km.shp","2 Kilometers" ,"FULL","ROUND","NONE"," #")

I need to have each buffered file to be named "WI_city_buffers*km.shp" with the * being the buffer size. How do I get the arcpy.Buffer_analysis to iterate through my list and apply each buffer distance and proper name?


texts: What are GIS methods for temporal analysis for poor data?

texts:
I'm interested to know what GIS methods or techniques there are for operationalising time-series data for spatial/temporal regression analysis where data have temporal data gaps.

For instance, if I am analysing survey data of a region, with survey data at 0, 5, 10 year snapshots, but want to employ some controls to assess temporal variations in data, what are some ways this can be done?


texts: From polygon to continuous measure of diversity?

texts:
First of all; I have tried to find a similar question, without success. Maybe it's because I'm quite new to GIS and I don't really know what exactly I'm looking for. If someone points me to a similar problem, I'd be happy to remove this post.

I need to create a 'continuous' or raster (in small grid cells) variable of population diversity for a given country. I have a shapefile showing the spread of the ethnic groups in polygons (fig. 1), and the outcome I'm looking for is a the 'mean indicator of diversity' in each of the administrative units (AU's, in this case, the 360 Nigerian constituencies).

Fig 1. Population groups polygons in Nigeria

The solution I came up with was to get the area percentage of each polygon in each AU, and calculate a heterogeneity index from that. But the problem is that I would be leaving aside quite a lot of information due to the distribution of administrative units. As shown in the fig. 2, the squares 'a', 'b', and 'c' would have the same 'segregation index', but it is clear that they are not in the same position vis-a-vis the 'hot spots'.

Fig 2.

So I though that another solutions could be create a grid map, and calculate the distance to the nearest border, but again sharing only one border is not the same as being in the central part of the map, where several groups live together.

After finding this question, I guess that the polygons could be transformed into points using their centroids, and then apply the same method. But the truth is that I'm new to this, and that question is not really clearly answered. How could I made such a thing?

Using another example, I want to create something like this (images from this website):

Given the distribution of some points with different qualitative characteristics, get a measure of diversity from where I could estimate a the 'mean heterogeneity' of each administrative unit.

How could I do it? I use R and QGIS, so I don't mind in which platform is the solution based on.


texts:
This past summer for research I was working on linking the NOAA Rivers Database (http://water.weather.gov/ahps/forecasts.php), The National Atlas Lake Database, and the the USBGN datasets so that there were shapefiles for each lake or river in the dataset instead of just one (or two) points for each river/lake.

We linked them in MATLAB and now have this done and plan to publish it open source. Our linking does a more complete job than what's currently in TNA for linking lakes and we think this could be useful to others in the community. However, we're not sure what's the best method to publish this so that the most people could use it. Do you think KML or geoJSON would be best? Or as they are currently, an array of links matching the data together?


texts: Schema Lock ERROR 000464 when ArcMap is open

texts:
I am having an issue with schema locks while working on feature classes within a file geodatabase in ArcMap.  There have been similar posts to mine (Avoiding exclusive schema lock error with ArcPy?).  However, I am not working in an ArcSDE environment.  I am running a Python script through the command line (Conda thru Windows if pertinent):
>>>python add_field_script.py
which calls arcpy functions WHILE ArcMap is open.
arcpy.DeleteField_management(fp_feature, ['field1', 'field2'])
I get this error (in command line console):
  File "field_mappings_config.py", line 42, in <module>
    arcpy.Delete_management(fc)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.6\arcpy\arcpy\management.py", line 4363, in Delete
    raise e
arcgisscripting.ExecuteError: ERROR 000464: Cannot get exclusive schema lock.  Either being edited or in use by another application.
Failed to execute (Delete).

My attrocious work-around is basically turn off and turn on again.  I CLOSE ArcMap and rerun.  Removing the fp_feature in the current map and/or collapsing feature dataset in ArCatalog without shutting off does not remove the Schema Lock.
Any ideas?


texts: Retrieve traditional Chinese entries from OpenStreetMap

texts:
I want to retrieve both simplified and traditional Chinese entries from OpenStreetMap. E.g.

<node id="244082315" lat="43.558611" lon="126.228611" version="5"     timestamp="2014-01-21T10:07:18Z" changeset="20119519" uid="274857" user="Supaplex">
<tag k="name" v="永吉县"/>
<tag k="place" v="suburb"/>
<tag k="gns:DSG" v="ADM3"/>
<tag k="gns:UFI" v="-1936843"/>
<tag k="gns:UNI" v="6771391"/>
<tag k="name:vi" v="Vĩnh Cát"/>
<tag k="name:zh" v="永吉县"/>
<tag k="gns:ADM1" v="05"/>
<tag k="name:zh_pinyin" v="Yongji Xian"/>
<tag k="name:zh-simplified" v="永吉县"/>
<tag k="name:zh-traditional" v="永吉縣"/>
</node>

Retrieves the simplified entry:

http://nominatim.openstreetmap.org/reverse?format=json&osm_type=N&osm_id=244082315&accept-language=zh

This does not work

http://nominatim.openstreetmap.org/reverse?format=json&osm_type=N&osm_id=244082315&accept-language=zh-TW

How can I retrieve traditional Chinese?


texts: QGIS, Number of Lines through a Point (of over lapping multiline shapefile)

texts:
I am doing an irrigation system analysis, for that there are multiple overlapping lines are passing through multiple points.

one number of line means a unit volume of the system,
and point means a gauging station.
and the number of lines passing means the flow throughput of the gauging station

How can I find the number of lines passing through a point?


texts:
I have been trying to find a marker animation plugin for OpenLayers to sync lat, lon, and time data along a track like this (in Leaflet):

http://leafletplayback.theoutpost.io/examples/example_2.html

Basically i want to draw a line string and animate the first marker till the tend with timeline slider with speed control.
I have seen : https://openlayers.org/en/latest/examples/feature-move-animation.html but it is confusing and doesn't meet my requirement.
Looking for a easy sample example.


texts:
I have a QGIS Project with many raster layers. I need to change/repair all paths since I moved the original data to Another disk.

layer = iface.activeLayer()

old_path = r'D:/Somefolder_QGIS_200317'
new_path = r'C:\GIS\Newfolder'

layer.setDataSource(layer.source().replace(old_path, new_path))

But im getting error:

  TypeError: QgsRasterLayer.setDataSource(): not enough arguments

What am I missing?


texts: Is it possible to change CRS without the .prj file?

texts: I'm working with QGIS and trying to change the CRS of some maps. Is is possible to do it without the .prj file?


texts: Load a feature from geoserver using WFS protocol on localhost

texts:
I am a student enginner in GIS , I want to load my features from geoserver using WFS protocol. I'm using this code but it's not appearing in my browser

    OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";
var map,wfs;

function init() {
//-11593508, 5509847, -11505759, 5557774
    var extent = new OpenLayers.Bounds(
        474813,  211897, 491713, 227436.5
    );
map = new OpenLayers.Map('map', {
    controls: [
        new OpenLayers.Control.PanZoom(),
        new OpenLayers.Control.Navigation(),
        new OpenLayers.Control.LayerSwitcher({'ascending':false}),
    ]
});
  var wms = new OpenLayers.Layer.WMS(
                    "test_postgres:territoire", "http://mylocalhost:8080/geoserver/test_postgres/wms",
                    {
                        LAYERS: 'test_postgres:territoire',
                        STYLES: '',
            format: 'image/svg',
                    },
                    {
                        isBaseLayer: true,
                    }
                );
  var wfsLayer = new OpenLayers.Layer.Vector("WFS", {
        strategies: [new OpenLayers.Strategy.BBOX()],
        styleMap: new OpenLayers.StyleMap({
            strokeWidth: 3,
            strokeColor: "#333333"
        })
    });

        map.addLayer(wfsLayer);

map.zoomToMaxExtent();


texts:
I need to save the symbology of several layers and I would like to do it with a PyQGIS-3 script.
I think I have to use the QgsMapLayerStyle class but I don't understand what is a QDocument and how to use this class
niv1 = QgsProject.instance().mapLayersByName('exempleNiveau1_UniteDiagnostic_CMH_valeur normalise')[0]
doc = QDomDocument("testdoc") # path to qml/xml document to create ?
elem = doc.createElement("maplayer") # layer with the symbology to save ? Vector layer is not supported
niv1.writeXml(elem, doc, QgsReadWriteContext())

doc = QDomDocument("C:/Users/Juliette/Documents/CRECQ/CartesIndicateurs/Layers/testdoc")
elem = doc.createElement("exempleNiveau1_UniteDiagnostic_CMH_valeur normalise")
niv1.writeXml(elem, doc, QgsReadWriteContext())


texts: WFS Stored Query is not stored correctly - it is stored with null in tag

texts:
If I want to list stored queries I get this error:

org.xml.sax.SAXParseException; lineNumber: 1; columnNumber: 735; The element type "null:ValueReference" must be terminated by the matching end-tag "</null:ValueReference>"

I found out, why this error occurs.

I created the Stored Query with this XML

<wfs:CreateStoredQuery service="WFS" version="2.0.0" xmlns="http://www.opengis.net/wfs/2.0" xmlns:wfs="http://www.opengis.net/wfs/2.0" xmlns:fes="http://www.opengis.org/fes/2.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:xsd="http://www.w3.org/2001/XMLSchema">
    <wfs:StoredQueryDefinition id="test">
        <wfs:Parameter name="gemarkungsname" type="xsd:string"/>
        <wfs:QueryExpressionText returnFeatureTypes="KRE_ALKIS:sg_flurstueck_f" language="urn:ogc:def:queryLanguage:OGC-WFS::WFS_QueryExpression" isPrivate="false">
            <wfs:Query typeNames="KRE_ALKIS:sg_flurstueck_f">
                <fes:Filter>
                    <fes:PropertyIsEqualTo>

    <fes:ValueReference>gemarkungsname</fes:ValueReference>${gemarkungsname}</fes:PropertyIsEqualTo>
                    </fes:Filter>
                </wfs:Query>
            </wfs:QueryExpressionText>
        </wfs:StoredQueryDefinition>
    </wfs:CreateStoredQuery>

But on the server in the data directory was this XLS stored. Here you can see <null:Filter instead of <fes:Filter.

<?xml version="1.0" encoding="UTF-8"?><wfs:StoredQueryDescription xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:fes="http://www.opengis.net/fes/2.0" xmlns:wfs="http://www.opengis.net/wfs/2.0" xmlns:gml="http://www.opengis.net/gml/3.2" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:xlink="http://www.w3.org/1999/xlink" id="test">
    <wfs:Parameter name="gemarkungsname" type="xs:string"/>
    <wfs:QueryExpressionText isPrivate="false" language="urn:ogc:def:queryLanguage:OGC-WFS::WFS_QueryExpression" returnFeatureTypes="KRE_ALKIS:sg_flurstueck_f">
        <wfs:Query wfs:typeNames="KRE_ALKIS:sg_flurstueck_f">
            <null:Filter xmlns:null="http://www.opengis.org/fes/2.0">
                <null:PropertyIsEqualTo>${gemarkungsname}<null:ValueReference>gemarkungsname</fes:ValueReference>
                </fes:PropertyIsEqualTo>
            </fes:Filter>
        </wfs:Query>
    </wfs:QueryExpressionText>
</wfs:StoredQueryDescription>

Can someone explain to me, why the null is saved instead of fes in my Stored Query in the data directory? I can edit the XML file. I can change here null to fes and everything works fine. But I like to understand why my Ajax not work correctly.


texts: Iterating/looping through Select by Attributes and Export Data?

texts:
In ArcMap 10.5.1, I have 2 feature classes in ModelBuilder; both have a common attribute "Well_Nu."

The goal is to loop for each record in the first classes's attribute table and do:

  1) Select by Attributes - Well_Nu = each value
  2) Export Data - Export: Selected features, with output feature class: each value_1mi.shp

And repeat this same process for the second feature class.

So I need to create individual feature classes for each (selected) feature.

How should this be done in ModelBuilder?

Noting that I would like to do it for the 2 feature classes, though it appears ModelBuilder only supports one iterator in a model; in this case, I could resort to creating 2 models.


texts: What does this goat shaped north-south line mean on Google Maps?

texts: If I zoom all the way out on Google Maps, just north-east of Australia is this goat-looking line that runs perpendicular across the equator. It's unlabeled.


texts: Automating relative filepath when adding new toolbox to arcmap catalog?

texts:
I'm creating an arcpy script for a custom toolbox and I want to make it so that a user can essentially utilize the script in that toolbox without having to change the source path (which I made relative in the properties of the script).

How do I go about doing so in Python or any other method?

I'm guessing it might be under the tool validation code somehow, but I'm not sure.


texts:
Everything I know about GIS I've learned myself.

I wrote a small utility for my non-profit employer that used the USGS 'The National Map's online services to pull a sized image of USGS Topographical Quad. I wrote the utility many years ago, and it was not only a crash-course in GIS, but in accessing a WMS.

The URL I had been using was similar to:

http://raster.nationalmap.gov/arcgis/services/Scanned_Maps/USGS_EROS_DRG_SCALE/ImageServer/WMSServer?service=WMS&request=GetMap&version=1.3.0&format=image/png&bbox=-105.6828,36.35569,-105.4633,36.45874&layers=0&width=1000&height=1000&crs=CRS:84&styles=default

..which would give an image of 1000px x 1000px around Taos, New Mexico. (Coordinates and image size would change, depending on user's area of curiosity.)

Now, however, that service was shut down. In my infinite glory, I had completely missed the posted notification and am now stuck without a running utility.

According to the USGS, the service has moved to 'services.arcgisonline.com/arcgis/rest/services/USA_Topo_Maps/MapServer'

So, learning what I could about the service, apparently WMS isn't available and I need to use the ExportMap functionality. This is the new URL I was able to translate the above to:

https://services.arcgisonline.com/arcgis/rest/services/USA_Topo_Maps/MapServer/export?bbox=-105.6828,36.35569,-105.4633,36.45874&bboxSR=EPSG:900913&layers=0&layerDefs=&size=1000,1000&imageSR=EPSG:900913&format=png&f=image

However, if you click on that, you'll see I get an awesomely helpful error of "bytes", code 500. (As a developer, I know Code 500 is Internal Server Error, of course..) But without a more helpful message, I'm not sure where I bungled it up.

Any thoughts?


texts: Altering the suggestQueryParameters in ESRI javascript API

texts:
I'm using a modified version of the CMV application and I'm trying to alter the query used by the suggestions query of a search widget. The suggestQueryParamters seems to allow for a custom where clause but I can't find any documentation on how to use it.

I want to alter the query from using wildcards on both sides of the search sting (the default) to only using it on the right side of the string i.e.

where field like '{value}%'
instead of

where field like '%{value}%'

In my code below you can see where I try to build a suggestQueryParams object into the sources property but it isn't having any effect. I'm guessing there is something wrong with the construction of the suggestQueryParams object but I can't find anything about constructing it correctly.

            this.searchDijit.on("load", function () {
            var sources = [
                {
                          featureLayer: new FeatureLayer("http://MyTestServer.amazonaws.com/arcgis/rest/services/clinrt/Building_XY/MapServer/0", {
                            outFields: ["*"],
                            //infoTemplate: new InfoTemplate("Parcels", "Number: ${ADDRESS}
Street: ${NAME}
Full_Address: ${Full_Address}")
                          }),
                          outFields: ["Full_Address"],
                          displayField: "Full_Address",
                          suggestionTemplate: "${Full_Address}",
                          name: "Parcels",
                          placeholder: "example: 1 CIVIC CENTER PLZ",
                          enableSuggestions: true,
                          maxSuggestions: 50,
                          suggestQueryParams: {where: 'upper(Full_Address) Like \'[value]%\'' }
                },

            ];

            //Set the sources above to the search widget
            this.set("sources", sources);
            this.set("activeSourceIndex", 0);
        });

Update
So while I couldn't find a way to alter the query that the search widget uses I did find a why to do what I wanted. I found that I could listen for the results event and alter the suggestion list box to remove the items I didn't want. It's not the prettiest or the most efficient way but it works:

            this.searchDijit.on( 'suggest-results', lang.hitch(
                this, function (e) {

                //console.log('suggest results', e);
                var itemi = this.searchDijit.suggestionsNode.children[0].children[0];
                if (!itemi.childNodes || itemi.childNodes.length == 0) return;
                //Here is the logic to remove unwanted items from the list of returned suggestions
                for (var i=0; li=itemi.childNodes[i]; i++) {
                    if(!isNaN(e.value) &&  li.innerText.indexOf(e.value) != 0){
                        itemi.removeChild(li);
                        i = i -1;
                    };
                }
            }
        ));


texts: Import shapefile to spatialiteDB and refer columns in QGis

texts:
Is there a way to import shapefiles into an existing spatialite-DB (or GeoPackage) and refer the columns of the importing shapefile to the right column of the DB.

E.g. I have a column in the DB named "abc". But in the shapefile the column is namend "def". Do I have to rename the columnname in the shapefile before import or is there a smarter way (especially I have a large number of shapefiles to import)?

Franz


texts: New GeoServer WFS error NoSuchAuthorityCodeException

texts:
What is causing this error?

org.opengis.referencing.NoSuchAuthorityCodeException: No code "EPSG:4326," from authority "European Petroleum Survey Group" found for object of type "IdentifiedObject".

https://geo.hyphae.design/geoserver/hyphaegeo/wfs?SERVICE=WFS&REQUEST=GetFeature&VERSION=2.0.0&TYPENAMES=hyphaegeo:PrescottGreenviewTest&TYPENAME=hyphaegeo:PrescottGreenviewTest&STARTINDEX=0&COUNT=10000&SRSNAME=urn:ogc:def:crs:EPSG::4326&BBOX=562073.76526905409991741,4185632.69297888176515698,563852.27126529649831355,4187445.67907886626198888,urn:ogc:def:crs:EPSG::4326:

Clicking on this link doesn't produce the error, but pasting it into the address bar does. More importantly, QGIS fails to load the layer with

2020-05-20T12:17:15 CRITICAL Layer Greenview test (Prescott) : Download of features for layer hyphaegeo:PrescottGreenviewTest failed or partially failed: Error when parsing GetFeature response : Error: not well-formed (invalid token) on line 1, column 1804. You may attempt reloading the layer with F5

I have verified:

the EPSG.data file is in its place
the WFS is on for the workspace
4326 is listed in the WFS settings page

Other than that settings are mostly OOB.


texts:
I have a Shapefile containing near 7000 pour points across the globe (GCS, WGS84). I'm using the following model to select the first feature point of the Shapefile, copy the selected feature in a new Shapefile, and then make its watershed. I need to do this because using the Watershed tool directly doesn't take account of the overlapping watersheds. Every thing is explained in my previous question.

The model works well, but it's painfully slow. The first two tools (Iterate and Copy) take less then a second to process, but the last tool (Watershed) takes up to 30 minutes for each watershed. In this question, I was advised to use Arcpy instead of Model Builder, but I've never used Arcpy or Python before.

So here's my question : Is there a site with templates for Arcpy where I can just fill in the blanks?

I'm looking at online tutorials, but since I'm in a hurry, having a pre-coded script would speed thinks up. All I need is a watershed script with a looping of the pour point used.


texts:
I am trying to use mapbox instead of leaflet to allow client to modify in the map if needed because I found in mapbox studio a friendly interface for modifying the map .

Is the free period in mapbox limited only for a month or is forever?


texts: Control histogram axis and legend with show_hist of rasterio

texts:
I have load raster (Sentinel 2) with Rasterio and I would like to create histogram of the pixels.
The problem is that I beleive the graph gives automatically the X label and also the legend labels.
This is how I create the histogram:
his_data=src.read([4,3,2])

from rasterio.plot import show_hist
plt.xlim(xmin=0, xmax = 1)
plt.xlabel('Reflectance')
show_hist(his_data, bins=200, lw=0.0, stacked=False, alpha=0.8,histtype='stepfilled', title="Histogram")

and this is what I get- the legend labels are not like the bands and the x axis called DN even though my image is reflectance and I defined it as reflectance:

is there any way to control this ? or I need to change it into array and then plot it


texts: How can I convert a polyline into a 10m wide polygon, like a thick polyline?

texts:
I want to convert a polyline into a polygon so that the polygon is like a thick version of the polyline.
E.g. if I had a centreline of a river and the river was 10m wide I would want to convert the polyline to a polygon covering the river. Is this possible in ArcGIS or QGIS?


texts: Viewing .grd and .gri file format in ArcGIS Desktop?

texts:
I have downloaded a file from DIVAGIS, the landcover map of India, and it is in the format of .grd and .gri file.

How can i view these files in ArcGIS 10.1?


texts: ArcGIS online - how to display a label show below 1:100

texts:
I'm trying to show cemetery label plots at all scales.  I have tried to change the visible range within the label features (see screenshot below) but this wont let me go below 1:100.  Anything that I type in below 100 will default back to 100.

The problem is that our cremation plots are approx 30cm x 30cm so even at a 1:100 scale these 3mm squares are too small to view labels (three digit numbers).

How do I allow labels to show below 1:100?
Is there an alternative way that would work?

EDIT:
I have allowed further zoom within the map settings of the application which allow for the map to zoom into 1:25 (yes the basemap will disappear at this level but that isn't a problem).  Hence why I should be able to see labels at this same scale.


texts: How to extract information from several shapefiles in the same time?

texts:
I'm using ArcGIS 10.2.1. I have a shapefile for 10 catchments and other three shapefiles for landuse, soil type and rock types. To extract the area of land use types, soil types or rock types for each catchment, I used Tabulate intersection.

However, Tabulate Intersection helps me to extract information from one shapefile only at a time (e.g. either soil type or land use type). I want to extract information from the three shapefiles in the same time.

I want to get the area of combination of land use, soil type and rock type in the same time.
The output I’d like to get should look like the table below

Any idea how can I do this in ArcMap?


texts: Removing two letters only if they appear at end of a field in QGIS labels

texts:
I am trying to remove the letters 'a' and 'o' only if they appear at the end of the field 'name' while labeling. For example, I'm working with names like:

Ontario, Manitoba, Alberta, Quebec, Nova Scotia

Effectively, I am trying to achieve something like:

Ontari, Manitob, Albert, Quebec, Nova Scoti

But if I were to use the expression
replace(name, 'a', '')
then all the 'a's would disappear like:

Ontrio, Mnitob, Albert, Quebec, Nov Scoti

Any ideas on how I should revise the code? I have already tried incorporating things like '%a' and right(name, 1) = 'a' to no avail.


texts:
I have coordinates and points from that coordinates. So what do I need to do to show point numbers on the map?

It's my second day with QGIS


texts:
I am new to ArcMap. I have an assignment to prepare map of each wards of a municipality. The data are point features. I need to print the map in different scales. like 1:2500, 1:5000 and also in different layouts like ISO A3, A1 landscape or portrait. The main aim of this assignment is to attain uniformity in symbology in different maps, to get standard font size and symbol size according to scale.

Questions:

Is there standards on symbol size, font size with respect to Paper Size?
I tried exporting layer (as .lyr file) and later imported in layer->properties->symbology. It works for me to make uniform symbols in different maps. Am I doing what I need to do or there is another better process?

Any reading material, link on these topic will be good for me.

Thanks,


texts:
So while using GeoServer 2.2 I did not have an issue with disappearing polygons in my display, now after upgrading to 2.7 my WMS now has scale dependent display of polygons, the data is there but not displaying till zoomed in. I suppose I have to make a new rule to add a SLD environment variable, maybe wms_scale_denominator? I am not overly comfortable with writing code yet and have looked at and tried to customize others examples with no avail. It should be simple but I am new and still learning.

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.1.0" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.1.0/StyledLayerDescriptor.xsd" xmlns:se="http://www.opengis.net/se">
  <NamedLayer>
    <se:Name>mainsorted</se:Name>
    <UserStyle>
      <se:Name>mainsorted</se:Name>
      <se:FeatureTypeStyle>
        <se:Rule>
          <se:Name>disturbance</se:Name>
          <se:Description>
            <se:Title>disturbance</se:Title>
          </se:Description>
          <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>dist_type</ogc:PropertyName>
              <ogc:Literal>disturbance</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <se:PolygonSymbolizer>
            <se:Fill>
              <se:SvgParameter name="fill">#6eb2ff</se:SvgParameter>
            </se:Fill>
            <se:Stroke>
              <se:SvgParameter name="stroke">#005ce6</se:SvgParameter>
              <se:SvgParameter name="stroke-width">1</se:SvgParameter>
            </se:Stroke>
          </se:PolygonSymbolizer>
        </se:Rule>
        <se:Rule>
          <se:Name>reclamation</se:Name>
          <se:Description>
            <se:Title>reclamation</se:Title>
          </se:Description>
          <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>dist_type</ogc:PropertyName>
              <ogc:Literal>reclamation</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <se:PolygonSymbolizer>
            <se:Fill>
              <se:SvgParameter name="fill">#55ff00</se:SvgParameter>
            </se:Fill>
            <se:Stroke>
              <se:SvgParameter name="stroke">#38a800</se:SvgParameter>
              <se:SvgParameter name="stroke-width">1</se:SvgParameter>
            </se:Stroke>
          </se:PolygonSymbolizer>
        </se:Rule>
      </se:FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts: Wrapping WMS across date line using GeoServer?

texts:
I am using GeoServer to serve a WMS version of the Natural Earth 1:10m Country shapefile.

I have everything working pretty well except that when it wraps around the dateline, there's an extra borderline in Russia and Antarctica:

I've read a few posts about how to deal with this in QGIS, but it seems like no matter what I do, there will always be a seam where GeoServer wraps the map.

Is there any way to get rid of these extra border lines like perhaps a way to selectively tell it not to draw a stroke between certain polygon vertices?

Any help would be much appreciated.


texts: Convert GeoTIFF height map to raw file for Unity terrain

texts:
I have a GeoTIFF file that QGIS has produced from merging multiple DEM TIFFs together. The gdalinfo for this file is:
Driver: GTiff/GeoTIFF
Files: ca_heightmap.tif
Size is 39601, 39601
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-125.000138888888998,43.000138888888898)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-125.0001389,  43.0001389) (125d 0' 0.50"W, 43d 0' 0.50"N)
Lower Left  (-125.0001389,  31.9998611) (125d 0' 0.50"W, 31d59'59.50"N)
Upper Right (-113.9998611,  43.0001389) (113d59'59.50"W, 43d 0' 0.50"N)
Lower Right (-113.9998611,  31.9998611) (113d59'59.50"W, 31d59'59.50"N)
Center      (-119.5000000,  37.5000000) (119d30' 0.00"W, 37d30' 0.00"N)
Band 1 Block=39601x1 Type=Float32, ColorInterp=Gray

How can I convert this to a raw file that I am able to import as a terrain in Unity? When I open this TIFF in GIMP, I just get a black/white image that loses the height map aspect (I'm not really sure what happens to it, but it looks like this in GIMP)


texts: In order to study land use and land cover change from LISS-III images I want to see the image coloured (i.e. false-colour composite RGB) using ENVI 5.0. What is the procedure?


texts: Find latitude of point on great circle given longitude

texts:
I am trying to find out wether my implementation has faults and/or if there is a better method.

The problem:

Given two points (Lat/Long) defined over the WGS84 ellipsoid, find the point lying on the great circle that connect them and has a given Longitude

My approach:

Convert the two points to ECEF coordinates
Derive the orientation of the plane where the great circle lies (cross product of the two ECEF vectors)
Derive the orientation of the meridian plane containing the point we are looking for (we know the longitude and that is perpendicular to the equator, so we only need sin/cos of the longitude)
Derive the ECEF vector of the sought point (is the intersection of the two planes, we need only the cross product of the two plane orientations)
Find the latitude of the point (is the arcsin of the third vector component, when normalized)

My doubts:

When converting to ECEF I was using the fact that the local earth radius is not constant, and thus I was computing it to correctly weight the three components. When I tried to plot my computed points vs the great circle on online maps (googlemaps, skyvector) they seemed a little bit off and I solved the discrepancy by removing the weighting (i.e., ECEF on a perfect sphere). Who is right?
Is there a better way? (that maybe avoids mapping to ECEF?)


texts: How can I read the settings for the QgsPaperItem from XML?

texts:
I saved a template qpt for the map composer and now I want to read all the settings. I am able to insert the items with composition.addItemsfromXML...
But this doesn't load the settings for the QgsPaperItem.
This is the code:

composition = newCompView.composition()
file = QFile(self.plugin.rootDirectory + "ui/PrintTemplates/" + self.printTemplate)
doc = QDomDocument()
doc.setContent(file, False)
docElem = doc.documentElement()
composition.clear()
paper = QgsPaperItem(composition)
paper.readXML(docElem, doc)
composition.addItem(paper)
composition.moveItemToBottom(paper)
composition.addItemsFromXML( docElem, doc)

With readXML the paperItem remains empty. I just managed to set it manually in the code. But I want this to be loaded from the template...


texts:
i managed to get my rgb-Image classified with Monteverdi. I discribed 8 classes with different classnames. For my work i dont need the rgb-values, i just need the classnames. Is there any way to export my classification result as vector file? Or can i save the information about the classname? I was just able to export my classification result as tif-file. So my input file contains 3-Band rgb values and my exported file has the same as well, just different values.

Iam open for all suggestions. Thanks.

Update: I am still not able to export my resulting classes. All i get is a 3-band rgb images regardless whether i choose *.img or *.tif. How can i do that?


texts:
i am trying to parse GeoJson from database with php but somehow i got an error and i cannot found any error in my script. here my script

$id = 1;
$db = new mysqli("localhost","root","","skriboo");
$sql = $db->query('SELECT * FROM profil_sekolah WHERE id_sekolah='.$id);
$sql2 = $db->query('SELECT * FROM pegawai WHERE id_pegawai='.$id);
$sql3 = $db->query('SELECT * FROM kelas_x A, kelas_xi B, kelas_xii C WHERE A.id_x=B.id_xi AND B.id_xi=C.id_xii AND C.id_xii='.$id);
$sql4 = $db->query('SELECT * FROM fasilitas WHERE id_sekolah='.$id);
$sql5 = $db->query('SELECT * FROM galeri where id_sekolah='.$id);
$sql6 = $db->query('SELECT bujur, lintang FROM lokasi WHERE id_lokasi='.$id);
$geojson = array(
'type' => 'FeatureCollection',
'features' => array(
    'type' => 'feature',
    'properties' => array(),
    'geometry' => array()
));

while ($row = $sql->fetch_array(MYSQLI_ASSOC)) {
    $properties = $row;
    array_push($geojson['properties'], $properties);
}

while ($row2 = $sql2->fetch_array(MYSQLI_ASSOC)) {
    $properties2 = $row2;
    array_push($geojson['properties'], $properties2);
}

while ($row3 = $sql3->fetch_array(MYSQLI_ASSOC)) {
    $properties3 = $row3;
    array_push($geojson['properties'], $properties3);
}

while ($row4 = $sql4->fetch_array(MYSQLI_ASSOC)) {
    $properties4 = $row4;
    array_push($geojson['properties'], $properties4);
}

while ($row5 = $sql5->fetch_array(MYSQLI_ASSOC)) {
    $properties5 = $row5;
    array_push($geojson['properties'], $properties5);
}

while ($row6 = $sql6->fetch_array(MYSQLI_ASSOC)) {
    $geom = array(
        'type' => 'point',
        'coordinates' => array(
            $row6['bujur'],
            $row6['lintang']
        )
    );
    array_push($geojson['geometry'], $geom);
}

header('content-type: application/json');
echo json_encode($geojson, JSON_NUMERIC_CHECK);
$db = null;

and then i check it in browser and i got this

array_push() expects parameter 1 to be array, null given in

maybe there is something wrong with my script and i hope there is an answer to solve this problem


texts: Matching vector network to slope raster dataset in ArcGIS Desktop Advanced

texts:
I am trying to calculate the average slope of a rivers network. Each part of the network is broken up into individually numbered segments. I have a high resolution slope raster and thought I could match this to the river network.

Is there an appropriate way to achieve this using ArcGIS Desktop?


texts: SQL filter more than 1 table in Cartodb.js

texts:
I have a button on my map where when user clicks, it will filter 1 cartodb layer with SQL. I want to allow each button click to filter two layers, each sharing the same field and value. Here is the working map

<link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/3.11/themes/css/cartodb.css" />

Supervisor Districts

All Districts

District 1

District 2

District 3

District 4

District 5


texts: Resolution changing between mosaic dataset and output from mosaic to new raster tool viewed in ArcScene?

texts:
ArcScene cannot accept the mosaic dataset (Raster 1) so I used the Mosaic to New Raster tool to create a raster (Raster 2) that we can use in ArcScene.

The problem which I face is the resolution in ArcScene is changed for Raster 2
but in ArcMap it's same as Raster 1.


texts: How to select features (by polygon) from multiple layers in QGIS?

texts:
I was trying to export features from multiple layers in a common region into new shapefiles after selecting them using the "Selection By Polygon" tool/button. However, the "Selection by Polygon" tool only selects features from the top/current layer (QGIS 3/2.99 as of July 2017).

Is there a way to select features from multiple layers (e.g. all visible layers or all layers) by drawing only one polygon?

I don't want to draw the query polygon for the region several times as it is tedious and the polygons I draw could be slightly different each time.


texts: Openlayers - Unable to set the style of features

texts:
I have tried to change the default style of my features, but nothing works. I have a database with the features of my user. These features I add manually to a Vectorlayer. But when i try to edit the style of those features. I have both tried feature.style = selected_polygon_style; and feature.setStyle(selected_polygon_style);. The first just doesn't change the style and the second one gives me a "AssertionError: Assertion failed." error.

let vectorSource = new VectorSource({
        format: new GeoJSON(),
        minScale: 15000000,
        strategy: bboxStrategy,
      });
      newFeatures.forEach((newFeature) => {
        let feature = new Feature({
          geometry: new Polygon([newFeature.geometry]),
        });
        feature.setId(newFeature.plotId);
        var selected_polygon_style = {
          strokeWidth: 5,
          strokeColor: "#ff0000",
          fillColor: "#aa2727",
        };
        feature.style = selected_polygon_style;
        vectorSource.addFeature(feature);
      });
      var vector = new Vector({
        //minZoom: 13,
        source: vectorSource,
      });
      this.setInteractionForPlotBoundriesLayer(vector, featureSelected);
      vector.set("name", "plotUserBoundriesLayer");
      this.plotsExtent = vectorSource.getExtent();
      this.map.addLayer(vector);

I use React combined with OpenLayers. So the this.map is created like as follows:

createNewMap() {
    this.map = this.createMap();
  }

  createMap() {
    return new Map({
      target: null,
      layers: [],
      view: new View({
        center: [594668.0262129545, 6602083.305674396],
        maxZoom: 19,
        zoom: 14,
      }),
    });
  }


texts: Programmatically creating a vector feature in OpenLayers2

texts:
I want to create a vector feature, programmatically, on my map. It's supposed to be in the A4-paperformat, for showing a sort of print preview.

I've created a control for handling this, but now I'm at the point where I need to create the feature. And I have no idea of how I can do it. Any ideas?


texts:
I have some troubles with GDAL module on Python. When I try to import gdal, console displays this:

Файл "c:/rstp/Lib/site-packages/osgeo/__init__.py", строка 48, из <module>    __version__ = _gdal.__version__ = _gdal.VersionInfo("RELEASE_NAME")
builtins.AttributeError: 'NoneType' object has no attribute 'VersionInfo'

or this with completely English translation:

File "c: / rstp / Lib / site-packages / osgeo / __ init__.py", line 48, from <module> __version__ = _gdal .__ version__ = _gdal.VersionInfo ("RELEASE_NAME")
builtins.AttributeError: 'NoneType' object has no attribute 'VersionInfo'

Then Wing moves to osgeo/__init.py file, where highlights this line:

__version__ = _gdal.__version__ = _gdal.VersionInfo("RELEASE_NAME")

Win10, Python version 3.7.4, GDAL 3.0.1

Can you help me please?


texts: Averaging XY points locations if in close proximity automatically over a large dataset in R

texts:
I have a point dataset (XYZ values) containing a few million points covering a 5km x 5km extent representing maxima locations. However there a great deal of these points which exist in 0.5m of others.

I would like to identify the points which are within 1m of one another and then average the coordinates to create a single point. I have written an R script to perform this task automatically on the XYZ data, see below. The process uses a loop, and operates on one point per iteration.

On large datasets this approach takes a long time to complete due to its iterative-loop nature, can any of you recommend a more efficient approach within R?

Thanks in advance,
Matt

##Sample data (17 point subset)
X <- c(671847.25, 671849.25, 671849.25, 671849.75, 671854.75, 671849.25, 671851.25,
671854.75, 671854.25, 671854.75, 671850.75, 671856.75, 671852.75, 671852.75,
671854.75, 671856.75, 671854.75)
Y <- c(3984125.75, 3984125.75, 3984125.25, 3984125.25, 3984123.75, 3984123.25,
3984123.25, 3984123.25, 3984121.75, 3984121.75, 3984120.25, 3984120.25, 3984119.75,
3984119.25, 3984118.75, 3984118.75, 3984117.25)
Z <- c(10.10000038, 12.31000042, 12.31000042,12.31000042,10.02000046, 11.30000019,
9.489999771, 10.02000046, 10.11999989, 10.11999989, 10.57999992, 10.11999989,
9.930000305, 9.930000305, 9.739999771, 10.72000027, 9.409999847)

maxima.data <- data.frame(X,Y,Z)

##set up empty output df to write new points to
results <- data.frame(x=as.numeric(), y=as.numeric(), z=as.numeric())
results <- rbind(results, c(0,0,0))
colnames(results) <- c("x", "y", "z")

while(nrow(maxima.data) > 0){
##Select first row in df, and perform a subset to identify all points within 1m
#always use first row
max.in.row <- subset(maxima.data, maxima.data[,1] < maxima.data[1,1]+1 & maxima.data[,1] > maxima.data[1,1]-1 &
maxima.data[,2] < maxima.data[1,2]+1 & maxima.data[,2] > maxima.data[1,2]-1)

##If more than 1 point/row returned, average the three columns XYZ to produce new
#point, otherwise the point stays the same and send to new 'results' DF
#And remove the point(s) in the original DF
#
if(nrow(max.in.row) > 1){
    #calculate column averages
    max.mean.row <- c(mean(max.in.row[,1]),mean(max.in.row[,2]),mean(max.in.row[,3]))
        #select rows in maxima.data to remove
    for(j in 1:nrow(max.in.row)){   #for loop no.1
        row.exclude <- which(maxima.data[,1] == max.in.row[j,1] & maxima.data[,2] == max.in.row[j,2])
        maxima.data <- maxima.data[-row.exclude,]
        rm(row.exclude)
        } #end for loop no.1

    #output
    results <- rbind(results, max.mean.row)
    rm(max.mean.row)

        } else {    #if only one row

    results <- rbind(results, maxima.data[1,])
    maxima.data <- maxima.data[-1,]

    }
    rm(max.in.row)

    }#end while loop

results <- results[-1,] #remove row of zeros
row.names(results) <- NULL


texts: Is there a quick way to invert a node selection in QGIS?

texts:
is there any way to invert a node selection in QGIS?

For example, I want to extract a small section of a long line - if I select the nodes I want to keep, is there any way to invert that selection so I can delete the unwanted nodes? The only way I can do it in v1.8 is to zoom all the way out, select all the nodes, zoom back in and deselect the ones I want to keep, then hit delete...which is a bit slow going...

Just wondered if there's a quicker way with a tool or keyboard shortcut?


texts: I am currently using rainfall data in determining rainfall patterns. But when I use Polygonize, the values on the attribute table of the created vector file were rounded off to the nearest whole number. I want to keep the decimal number but iI can't find any option to do that.


texts: How to detect when new markers have been added to the leaflet map?

texts:
May i ask how to I add a event listener to the map to detect when a marker has been added to the map?
In the docs, it says that we can use the "add" event, but how do I actually use it?
http://leafletjs.com/reference.html#marker

I have tried:

map.on(marker,'add', function() { // Throws a error "marker is not defined"
    console.log('marker added')
});

map.on('add', function() { // doesn't do anything
    console.log('marker added')
});

Any help?


texts: How to compare points joined to grid, versus points joined to grid + 15km buffer, using sf package?

texts:
I have an issue where I'm spatially joining points to grid squares, but sometimes they're just barely out of range. To address this, I want to identify which points become incorporated into these original grid cells, when applying a buffer of 15km.

What is the best way to go about this using the sf package? How do I create a new buffered-square join I can directly compare to the original squares? Particularly since this new buffered join will overlap several other grid cells, such that many points will be double or quadruple counted?

At the end of the day, I'd like to examine an object or dataframe that identifies every original square and the points that are 'missed' just outside it, within that 15km border.

library(tidyverse)
library(sf)
gridpolygons <- read_rds('Data/gridpolys.rds') %>% st_as_sf
gridpolygons.metersproj <- st_transform(gridpolygons, 32119)

# create 15 km buffer and transform back to wgs84
buffer15km <- st_buffer(
  gridpolygons.metersproj, dist = 15000) %>% st_transform(4326)

points <- st_as_sf(getPoints()) # Pulls shapefile of points into sf object
                                # (custom fn; then conversion to sf object)

# Original join
cell_points <- st_join(gridpolygons, points)

# Buffer object join? e.g.,
  # buffer_points <- st_join(buffer15km, points)#?
# Or alternate method to avoid double-counting, that I can compare
  # compare easily to first spatial join?
# Method that cuts a hole out of the buffered polygon
  # to make a polygon that's just border, and then joins points with replacement?

The objective is to identify the points that are 'near misses' from being joined to each cell, or within 15km of the grid cell.


texts: Hard line-barriers for Thiessen polygons

texts:
Similar to Creating polygons from points using lines as barriers, which has not been resolved; plus, the scenario is slightly different.
I have one dataset with 60,000 points and one with boundaries of ca. 5000 districts.
I need to establish areas for the points (using Thiessen).
However, the resulting areas must not cross district boundaries.
As a first step, I'd just created the Thiessen polygons. Then I'd make a union with the districts. I could also delete (or dissolve) the residuals that contain no points. (in the example image: black lines)
But: this actually solves nothing. In the example image, the residual of the western district (encircled in blue) would have to be divided up between points A and B. I can't think of any solution to this.


texts: How to clip polygons with a fishnet layer?

texts:
I have a layer with multiple polygons (input) and a fishnet layer in polygon (clip layer).

I want to clip the multiple polygons with the fishnet to be able to calculate how big portion of the input layer that cover each cell in the fishnet layer.

If I select a fishnet cell and run the clip tool I get what I want, the area of the input polygon(s) covering the fishnet cell.
But if I dont have any cells selected from the fishnet layer and want to run it with all of them, it does not work, it clip the layer with the extent instead of the independent polygons in the fishnet.

Any suggestions of how to do this? I have arcgis 10


texts: How to use only English in names with a local OpenStreetMap server?

texts:
I have a local OpenStreetMap Server (Ubuntu + postgresql + mod_tile + mapnik).
How to use only English language in names (or any Specific language) for all areas?

P.S. I understand that names store in tags and maybe for some areas it's not available.

UPDATE Can I change style file for mapnik to render tiles in specific language?


texts: QGIS/ArcGIS Pro: WMTS-Services - different display quality

texts:
Using ArcGIS Pro (1.4) and QGIS 2.18 I discovered that both programs have a different "display-quality" of the same WMTS-Service (in my case basemap.at - the map called basemap-grey). I was surprised and checked some properties, but with no success. ArcGIS Pro seems to render it "smoother" (Screenshots attached). In both programs I tested it with the same extent at the same scale (and EPSG 3857 as data-frame SRS). In QGIS I manipulated some resampling-parameters, but without getting the same quality.
Does someone know how to improve it ?


texts: How to remove clouds in High resolution Imagery(.5 mGeoeye)

texts:
I have imagery with high resolution .5m, large parts covered with clouds, so it's obscure many urban areas.

How can I remove the clouds?


texts: What needs to be done to downgrade a ESRI addin from 10.1 to 10?

texts:
I have created an esri addin which targeted ArcGIS 10.1, however I need to have the addin to work on ArcGIS 10, what needs to be done within my solution to get it working?
I have checked the references, I have changed the targeted version witihin the config.xml, but when I try and run the addin (from Debug mode in VS2010) I don't see the addin as I would on the 10.1 version.


texts: What is the canonical way to perform raster calculations using R stars?

texts:
I have a data cube that I am handling in the R stars package. It has 3 dimensions (latitude, longitude and time) and 2 attributes (northwards and eastwards components of wind). A dummy equivalent dataset can be generated as follows:
library(units); library(stars)

(w <- system.file("nc/bcsd_obs_1999.nc", package = "stars") %>%
  read_stars() %>% drop_units() %>% setNames(., c("u10", "v10")))

I want to calculate the resultant magnitude of wind for every pixel using simple trigonometry. My data is considerably larger (87m pixels), so efficiency of computation is very important
What is the canonical way to perform such raster calculations?
Ideally I'm looking for answers that are generalisable to any raster calculation problem using R stars, not just the specific example I give here - mostly because I know I'll need to be doing more complicated calculations in the future.
I've found three options:
Method X
system.time(
  x <- sqrt(w$u10^2 + w$v10^2)
)
# 0.001s
class(x)
# [1] array
x <- st_as_stars(x, dimensions = attr(w, "dimensions")) # Hack it back to a stars object

This is the most intuitive method (at least to me) and is very fast - but the resulting object has lost the dimensions and is no longer a stars object, so presumably this isn't the canonical method.
Method Y
system.time(
  y <- sqrt(w[1,,,]^2 + w[2,,,]^2)
)
# 0.006s

This generates a stars object but takes 6 times as long as the previous method and is less intuitive.
Method Z
wm <- merge(w)
system.time(
  z <- st_apply(wm, 1:3, function(ws) sqrt(sum(ws^2)) )
)
# 0.167s

This seems the most canonically "stars" method, but requires reorganisation of the data and fitting a function that could get complicated and prone to errors if computing something involving many attribute layers. Plus it's super slow even discounting the required merge operation to convert the attributes to a dimension.


texts: Geocoding through Canadian French language barrier?

texts:
This is related to Very Inaccurate Geocodes in Canada (Quebec).

It has come to my attention that the streets file I use as my address locator in Canada uses French (or simply accented) characters for the names. These show up in similar manners to the following:

QuÃ©bec
De la RiviÃ re-du-Lac
FranÃºois-de Villars

The Address product, however, uses standard (Oops, my Ameri-centrism is showing?) English characters, causing the previous examples to display like so:

Quebec
De la Riviere-du-Lac
Francois-de Villars

I believe the geocoding would fair much better if the street names actually matched via the same character set. How would you recommend I go about converting from one to the other? I have considered regex'ing it, but I thought I would ask around before I start.

P.S. It would be nice if I didn't get garbage text from French characters, regardless. All the info is stored in Geodatabases, by the way, if that matters at all. I'm on Arc10, as well.

P.P.S. It would appear that I can display accented characters just fine, but in the Shapefiles and GDB Feature Classes, they become garbage text. Is there some setting I am missing here?    CURIOUS: It would seem that in one street product the names come up garbage, but in the other they display just fine.


texts: Optimize OpenLayers2/ Geoserver 2.5 with GeoWebCache - reprojection necessary?

texts:
I have an OpenLayers map setup with several base layer options (some served by my own GeoServer) and many optional layers served by my own GeoServer.  I've selected Google's EPSG:900913 as the map's projection.  And I've specified specific resolutions to optimize use of the GeoWebCache (see below for explanation).  Here are the relevant portions of the map's options when created in OpenLayers:

 options = {
    projection: new OpenLayers.Projection("EPSG:900913"),
    displayProjection: new OpenLayers.Projection("EPSG:4326"),
    numZoomLevels: 12,
    resolutions: [152.8740566, 76.43702829, 38.21851414, 19.10925707, 9.554628536, 4.777314268, 2.388657134, 1.194328567, 0.597164283, 0.29858214, 0.14929107, 0.07564554],
 ...

These resolutions come from looking at the Tile Matrix Set for GridSet EPSG:900913 in GeoServer.  These correspond to zoom levels 10 through 22.  I did this to optimize the tile requests to match the seeded cached tiles in GeoWebCache.

Turns out that many of my layers, served by my GeoServer, have native projection of EPSG:2234 (NAD83 / Connecticut (ftUS)).  I see that I can, if I wish, create a new GridSet for caching tiles in this projection.  If I just let it create default Tile Matrix resolutions, they will be different from the resolutions for EPSG:900913 GridSet.  At least I believe so - these will be in feet rather than in meters, so I'm comparing apples to oranges anyway. But now I'm thinking this whole step is completely unnecessary...

Can I still simply seed my native EPSG:2234 layers with the EPSG:900913 Tile GridSet, and then expect tiles requested by OpenLayers to the GWC service using srs:'EPSG:900913' will get created and cached properly?  Or should I force reprojection at the layer level, so it is both published and cached at the same projections?  I don't want to go through this step if it is completely unnecessary...


texts: Failing to Show Custom Sentinel 2 Image Collection on a Map in Google Earth Engine

texts:
I'm creating a custom image collection over a bound region as follows:

Get image collection containing images taken by satellite in 7 days

var collection_to_covert_to_image = ee.ImageCollection("COPERNICUS/S2").filterDates(d1, d2).filterBounds(roi);

Convert image collection to a single image:

var img = collection_to_convert_to_image.toBands();

Now this creates an image with bands list as a sum of all images in
collection. Like instead of correct B1, B2, B3, B4, etc it is creating
bands list having IMAGE_1_DATE_B1, ...IMAGE_1_DATE_B2, ...,
IMAGE_2_DATE_B1, IMAGE_2_DATE_B2, ...
What I want the resulting combined image from image collection to have
B1, B2, B3, B4, ...

Repeat Steps 1 and 2 to get multiple 'combined' images and then creating image collection as follows:

var custom_col = ee.ImageCollection([ combined_img_1, combined_img_2, ... ]);

Finally to preview this collection, I add first image to the map:
Map.addLayer(custom_col.first(), { min: 0, max: 6000 });

But nothing gets displayed although layer is loaded. I believe this is because of new images not having proper B1, B2, etc bands
So how to make combined image that contains these bands as sum of child images' respective bands like e.g combined_img_B1 = img1_B1 + img2_B1 + ... in order to display it on map?


texts: Does anyone display OpenStreetMap through MapServer, ArcGIS for Server or Geoserver?

texts:
I use Windows OS. The data (OpenStreetMap) are stored in a postgres/postgis database.

How do I display this data correctly?

Update: Now I use: geoserver-2.3.5.exe, postgresql-9.1.9-1-windows-x64.exe, postgis-pg91x64-setup-2.0.3-2.exe, osm2pgsql

Update2 I lookd at the source for "osminabox". https://git.hsr.ch/index.html#repositoryBrowser;26067c88-f4a8-46e8-be0d-95c1a5091c94;null;Code/osm2gis/geoserver/geoserver_data/styles

But I find many sld styles (I have only 4 layers: planet_osm_polygon, planet_osm_line, planet_osm_point, planet_osm_roads).
There are a lot of styles for points: osm_poi.sld, osm_building.sld, osm_pow.sld. I don't have this layer.
Also I find strings such: "img/recycling.png"

So, If I use intranet system I can't work?

Update 2 continue:
I try to use https://github.com/fotographs/OSM-Styles-for-GeoServer/tree/master/sld. But I cannot find sld styles for the points.

https://help.openstreetmap.org/questions/7521/geoserver-style-files - Would it help?


texts: GDAL retile to create near-uniform file sizes?

texts:
I have a 60Gb raster vegetation dataset. I need to be able to access single pixels relatively quickly. At the moment, I am retiling at a fixed size, and saving those, indexing them in Python, and accessing them from there.

However, this results in tiles with a wide array of sizes - places with consistent vegetation get compressed well, and are really small, places with a lot of variance are much bigger. I would ideally like to optimise the file sizes so that my index is not too huge (~10K files maybe), but so that file access is also minimised - for any given pixel, I would rather not try and load a 100Mb file.

Is there an easy method of retiling in one step that would tile a dataset into (very) approximately uniform file sizes? I guess I could tile at a large size, and then re-tile each of the largest files into smaller files, and repeat. But it would be nice if there was a way that didn't require multiple steps.


texts: OpenLayers Style Question ... not working in IE

texts:
I'm trying to use OpenLayers to produce a map I'm using the following to style a layer in the map.

It seems to work fine in Firefox but I get javascript error in IE 'symbolizer' is null or not an object ....?
Any Ideas?

var previewVectors = new OpenLayers.Layer.Vector("previewVector", { isBaseLayer: true,

          styleMap: new OpenLayers.StyleMap({
            "default": new OpenLayers.Style({
              strokeColor: "white",
              strokeWidth: 1
            }, {
              rules: [
                new OpenLayers.Rule({
                  filter: new OpenLayers.Filter.Comparison({
                    type: OpenLayers.Filter.Comparison.EQUAL_TO,
                    property: "leUsage",
                    value: "Bulk Commodity"
                  }),
                  symbolizer: {
                    fillColor: "olive"
                  }
                }),
                new OpenLayers.Rule({
                  filter: new OpenLayers.Filter.Comparison({
                    type: OpenLayers.Filter.Comparison.EQUAL_TO,
                    property: "leUsage",
                    value: "Service Industry"
                  }),
                  symbolizer: {
                    fillColor: "pink"
                  }
                }),
                new OpenLayers.Rule({
                  filter: new OpenLayers.Filter.Comparison({
                    type: OpenLayers.Filter.Comparison.EQUAL_TO,
                    property: "leUsage",
                    value: "Car Precincts"
                  }),
                  symbolizer: {
                    fillColor: "blue"
                  }
                }),
                new OpenLayers.Rule({
                  filter: new OpenLayers.Filter.Comparison({
                    type: OpenLayers.Filter.Comparison.EQUAL_TO,
                    property: "leUsage",
                    value: "Container Trade"
                  }),
                  symbolizer: {
                    fillColor: "green"
                  }
                }),
                new OpenLayers.Rule({
                  filter: new OpenLayers.Filter.Comparison({
                    type: OpenLayers.Filter.Comparison.EQUAL_TO,
                    property: "leUsage",
                    value: "Warehousing and Container Park"
                  }),
                  symbolizer: {
                    fillColor: "navyblue"
                  }
                }),
                new OpenLayers.Rule({
                  filter: new OpenLayers.Filter.Comparison({
                    type: OpenLayers.Filter.Comparison.EQUAL_TO,
                    property: "leUsage",
                    value: "Other"
                  }),
                  symbolizer: {
                    fillColor: "brown"
                  }
                }),
            ]
            })
          })
        });


texts: How to get mean geographic direction of a polygon in QGIS?

texts: I have a raster height model and would like to calculate the horizontal geographic direction (north, east, etc. - orientation) of each cell, then calculate then the mean value of all these cells in one region (polygon). This is to get information about to what direction a vineyard is exposed to. How do I do this?


texts: SAGA Raster Mosaic Issue with output (projection/not showing)

texts:
I am trying to create a mosaic from a couple of georeferenced raster images (old aerial photos) using the SAGA mosaic raster layers function from the QGIS toolbox. Project and layers are all on the same CRS (EPSG:4674). I kept the default settings from the SAGA Mosaic Raster Layers window and processing runs smoothly. However, the mosaic output layer (named 'Grid') does not show on the project CRS nor when I "go to layer". In properties, the output's CRS is the same as other layers and project.

Any ideas of what might be happening, and how to fix it?


texts: Specifying EPSG:31370 CRS in R Leaflet for WMS-service returns blank tiles

texts:
I am trying to load historical map tiles from a WMS-server (GetCapabilities) and display them using R Leaflet.

Using the default CRS in Leaflet works fine, using e.g.:

leaflet() %>%
  setView(4.69839084, 50.880089, zoom = 14) %>%
  addWMSTiles(
    "https://geoservices.informatievlaanderen.be/raadpleegdiensten/histcart/wms",
    layers = list("ferraris"),
    options = WMSTileOptions(format = "image/png"),
    attribution = "Source: GIS Geoservices Informatie Vlaanderen")

However, the projection should be Belge 1972 / Belgian Lambert 72 (EPSG:31370).

I can specify an example CRS (EPSG:3006), and the tiles still load:

epsg3006 <- leafletCRS(crsClass = "L.Proj.CRS", code = "EPSG:3006",
  proj4def = "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
  resolutions = 2^(13:-1), # 8192 down to 0.5
  origin = c(0, 0)
)
leaflet(options = leafletOptions(crs = epsg3006)) %>%
...

However, using the epsg.io Proj4js-string for EPSG:31370 results in gray/blank tiles in R Leaflet:

epsg31370 <- leafletCRS(
  crsClass = "L.Proj.CRS", code = "ESRI:31370",
  proj4def = "+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs",
  resolutions = 2^(13:-1), # 8192 down to 0.5
  origin = c(0, 0))

leaflet(options = leafletOptions(crs = epsg31370)) %>%
  setView(4.69839084, 50.880089, zoom = 14) %>% # center on Leuven
  addWMSTiles(
    baseUrl = "https://geoservices.informatievlaanderen.be/raadpleegdiensten/histcart/wms",
    layers = list("ferraris"), # request Ferraris-tiles from available WMS tilesets
    options = WMSTileOptions(format = "image/png"),
    attribution = "Source: GIS Geoservices Informatie Vlaanderen") %>%
  addControl("Leuven on the Ferraris map (1777)", position = "topright")

Additionally/possibly related:

I'm not sure how the resolutions-arguments work. Should I be able to find that through WSM GetCapabilities?
Proj4js returns a depreciation warning:

  Warning in leafletCRS(crsClass = "L.Proj.CRS", code = "ESRI:102003",
  proj4def = "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0
  +y_0=0 +datum=WGS84 +units=m +no_defs", : projectedBounds argument is
  deprecated and has no effect, use the bounds argument.


texts:
I activate the TILED parameter to True in a WMS Layer and I get the following strange behavior as shown in the images:

This happens only for a couple zoom levels. If I zoom out/in more the features are rendered in the correct position (as in the first image).

Any idea why is this happening?

Here is my code for the layers:

        // layer komvoi
        wms_layer_komvoi = new OpenLayers.Layer.WMS( "Όλοι Κόμβοι","http://localhost:8080/geoserver/wms",
            {layers: 'rural:komvoi_real',
            transparent: true,
            TILED: true
            },
            {
            format: "image/png",
            visibility: false,
            displayOutsideMaxExtent: false,
            projection: new OpenLayers.Projection("EPSG:900913")});
        // layer network
        wms_layer_network = new OpenLayers.Layer.WMS( "Ολόκληρο Δίκτυο", "http://localhost:8080/geoserver/wms",
            {layers: 'rural:foc_network',
            transparent: true,
            TILED: true
            },
            {
            format: "image/png",
            visibility: false,
            displayOutsideMaxExtent: false,
            projection: new OpenLayers.Projection("EPSG:900913")},
            {
            });

        // add layers
        map.addLayers([wms_layer_komvoi,wms_layer_network,wms]); // first you put the layer with the points
        map.zoomToExtent( new OpenLayers.Bounds(36.190737, 35.353912,35.150577, 50.274810));
        map.setCenter(new OpenLayers.LonLat(22.228697,37.650344), 8);


texts:
I have a month's worth of snowfall contoured polygons in one feature class. Polygons represent one day and has a field called SNOWFALL. How can I display the total snowfall for the month? I have the standard license.

I have attached a screenshot of the identity tool on a polygon.

I tried a union, then after that a dissolve by Shape area with Sum of Snowfall checked. However many polygons have the same area so that didn't work.


texts: Getting watermarks on my satellite image

texts:
I am Using QGIS 2.16.1. When I open the Google satellite image by using open layers plugin I am getting watermarks on it. So many time I uninstalled and installed. Even I installed new version of QGIS. I went through the same question replied by you.
How to close watermarks window in QGIS Openlayers plugin?

but the link you suggested has been removed already.