texts:
I am trying to document a code for a QGIS plugin. I have installed sphinx in anaconda. When I try do make html, i get an error saying "No module named 'qgis'", I belive it comes from one of the classes trying to import qgis.core. I have read suggestions to use OSGeo4W to install QGIS as it would give me all the enviroment there. However, when i try to import qgis.core within that shell it still cant fin qgis module. I have also faild to install sphinx in OSGeo through pip.

So my question is, how can I use sphinx to document a QGIS plugin code?

I am using anaconda 5.1, python 3.6.5 and sphinx 1.7.2


texts: Joining table to feature class AND matching fields using ArcMap

texts:
I am digitizing points from pdfs and aerial imagery.  Just points.  Then I want to create lots of fields and populate them without resorting to FieldCalculator or arcpy.DA modules.  There are too many fields and rows to add manually using FieldCalculator.  So here are my first steps:

create new feature class --> point
manually add (digitize) 30-40 points

Then I need to populate all 30 rows with 10 fields worth of values. So I am separately filling out a csv with 30 rows to join to those points.  For instance 30 features/points on new feature class gets 30 rows in a csv. The csv has 10 attribute fields and I fill out all the values for each of the 30 rows ON the CSV because it's fast and easy.  Next steps:
3. join the csv to the basically empty feature class
4. export the joined feature class
Then each week I digitize 40 new points in feature class and I need continue with the csv protocol and join.  So I tried:
5. edit existing feature and digitize the 30 new points
6. create new csv with id_join field matching the objectid of 30 new points in feature class (31-60 for example).
My problem is, apparently there is no capability to essentially join WHILE mapping the attributes from csv to the feature class.  For instance if field: num_people was created on first join and existed in feature class when second join occurs, I can't join my batch with rows 31-60 containing the EXACT same fields and field names and populate those fields from 31-60.  Brand new fields are created.
So it's basically a join/append hybrid I am trying to achieve to allow a periodic update of growing feature class.
Does anybody have a workflow preferably using joins or other tools in ArcMap?


texts: Grouping country polygons into a region without dissolving borders

texts:
PowerBI allows for the import of a shapefile via JSON to create a map for data visualization.
In PowerBI, I want my users to be able to select either a single country or a region that contains that country.
As far as I can see, my process in QGIS should be:

Create a shapefile with all country polygons
Create a shapefile in which I group the country polygons into regions and delete polygons outside the regions
Merge the two shapefiles
Convert to JSON

If this process is correct (please tell me if there is a better way!), my problem is that using Merge Selected Features deletes the borders between the countries in the region, resulting in a very politically incorrect map!

Is there a way to group country polygons into a region without deleting borders?


texts: Problem refreshing GroupLayers in ArcScene (ArcObjects 10 .NET SDK)

texts:
The functionality I am seeking is to add a FeatureLayer (containing some data from a file selected for import) inside of a GroupLayer collecting all of the data from that particular source. The logic for creating the feature layer and adding it to the scene, and subsequently forcing it to refresh so that it draws, has been working for FeatureLayers some time. However, it doesn't work for GroupLayers. The refresh method is as such:

    public void RefreshLayerInDocument(ILayer layer)
    {
                var sxdoc = this.Application.Document as ISxDocument;
                sxdoc.Scene.SceneGraph.Invalidate(layer, true, true);
                sxdoc.Scene.SceneGraph.Invalidate(sxdoc.Scene.SceneGraph.ActiveViewer, true, false);
                sxdoc.Scene.SceneGraph.RefreshViewers();
                sxdoc.UpdateContents();
    }

Note that manually selecting Refresh on the TOC allows a proper refresh and everything works.

The code for creating the grouping and adding the feature layer to it takes the form;

            IGroupLayer newGroupLayer = new GroupLayer();
            newGroupLayer.Name = layerName;
            scenehook.Scene.AddLayer(newGroupLayer);
            newGroupLayer.Add(layer);
            RefreshLayerInDocument(group as ILayer);
            RefreshLayerInDocument(layer);

Altering the order of the Refresh calls on the group and the featurelayer has had no impact on the outcome, it still requires a manual refresh.

My thoughts are that there must be an invalidate call required for groups that I'm not aware of, but I haven't had any luck searching the documentation.


texts:
When I run QGIS the log messages panel shows something like this:
problem with SAGA installation: SAGA version non supported (found 7.4.0 required 3.2.). But in a ubuntu shell the command "saga_cmd" returns SAGA Version 2.3.1.
I removed and reinstalled SAGA but the problem is the same and isn't possible to use SAGA. Any idea?


texts: Working with irregularly spaced gridded NetCDF data in R

texts:
I'll be working with modelled air pollution data soon. I've been given some sample data and it looks like the grid is very irregular ('brick' will throw an error about irregular grid), which makes it harder to get it into raster format.
How do I get this irregular data into a spatial object? (SpatialPointsDataFrame or raster)
This is what I tried so far:
# load packages ####
library(ncdf4)
library(raster)
library(rgdal)
library(fields)

# netcdf to data matrix ####
nc_data <- nc_open('filename.nc')

cnc_PM2_5       <- ncvar_get(nc_data, attributes(nc_data$var)$names[1])
longitude       <- ncvar_get(nc_data, "longitude")
latitude        <- ncvar_get(nc_data, "latitude")
#time           <- ncvar_get(nc_data, "time")

# plot matrix, latitude is upside down ####
image.plot(longitude,rev(latitude), PM2_5,
           main="PM2_5", ylab="latitude")

# get quick shapefile for Finland (http://www.diva-gis.org/datadown)####
finland <- readOGR("FIN_adm4.shp")
# add map projection
proj4string(finland) <- CRS("+proj=longlat +datum=WGS84 +init=EPSG:4326")
helsinki <- finland[finland$NAME_4=="Helsinki",]
# add Helsinki outline to image.plot
plot(helsinki, add=TRUE)

Reproducible example
    PM25 <- matrix(data=c(4.540674,3.437939,3.373072,2.335204,2.140603,2.529804,
                          4.475807,4.346074,6.42181,3.113605,2.594671,2.854138,
                          3.632539,6.097476,3.308205,5.643409,5.643409,2.140603,
                          4.151473,5.643409,2.529804,2.400071,3.827139,4.994741,
                          3.243339,3.373072,2.075737,2.400071,4.281207,3.697406,
                          2.854138,5.448809,2.075737,3.827139,3.567672,2.854138,
                          4.475807,2.983871,3.697406,5.578542,2.140603,2.20547,
                          6.22721,3.502806,5.773142,4.994741,5.254208,2.594671,
                          3.502806,3.827139,3.502806,2.400071,3.762273,2.724404,
                          3.892006,3.502806,2.724404,2.20547,4.151473,5.838009,
                          6.746144,5.059608,2.140603,2.270337,3.567672,2.140603,
                          3.437939,2.724404,2.140603,3.437939,3.308205,4.994741,
                          4.086606,2.335204,2.140603,5.578542,4.994741,2.724404,
                          3.892006,4.605541,4.281207,4.281207,2.140603,2.529804,
                          3.308205,3.892006,3.697406,3.308205,2.270337,6.097476,
                          3.308205,4.41094,5.319075,4.086606,2.854138,6.162343,
                          2.140603,3.308205,3.048738,4.02174), ncol=10, nrow=10)

    # latitude
    latitude <- c(60.1191177368164,60.1192321777344,60.1193504333496,
                          60.1194686889648,60.1195831298828,60.119701385498,
                          60.1198196411133,60.1199340820312,60.1200523376465,
                          60.1201667785645)
    # longitude
    longitude <- c(24.5949993133545,24.595235824585,24.5954704284668,
                          24.5957050323486,24.5959415435791,24.5961761474609,
                          24.5964126586914,24.5966472625732,24.5968818664551,
                          24.5971183776855)

    image.plot(longitude,latitude, PM25, main="PM2_5")


texts: Shapely's simplify won't simplify any more


texts: How to reset cql_filter of WMS on OpenLayers 3?

texts:
I can dynamically change the CQL filter, but I want to reset it to default (get WMS without any filter).

Class I am using  "ol.source.TileWMS" as source from GeoServer

Here is my code to which I used to update the CQL filter:

overlay_layer_definition_array[selected_layer].getSource().updateParams({  'LAYERS': Lay_Param[selected_layer][1] + ':'
        + Lay_Param[selected_layer][2]   ,'cql_filter': ''+strUser+strUser1+strUser2});

The above code works fine, but now I don't want to pass any filter, tried all of below

How to immediately refresh ol.source.TileWMS  does not work

refresh() does not work

'cql_filter': nil does not work

'cql_filter': 'nil'  does not work


texts: Extracting raster DEM elevations to vector attribute tables in ArcGIS Desktop?

texts:
I have an ArcMap 10 project with multiple layer groups. One layer is a raster DEM, another layer group is a bunch of vector shapefiles (including points, lines and polygons). The client wants all vector data to contain Z values.

Populating the attribute tables for the vector data with X and Y values is easy enough but how do I populate the Z values?

Is there a way that I can use Field Calculator to "grab" the height values from the DEM layer for each vector feature?

I should add that we don't have a licence for the 3D Analyst or Spatial Analyst extensions.


texts: Trying to perform kernel density analysis but obtain: 'error all input points identical' or 'error in executing grid expression'

texts:
I am trying to perform a density analysis of a point shapefile on ArcMap. I am unsure of what to insert for population field. I've tried the actual number of people in the census tract where my point is and also tried adding a column to insert count (of points), so that each point has a count of 1. I have also left it as None. Neither works.

I get 2 types of errors:

all input points are identical
errors in executing grid expression


texts: Merging shapefiles with the same name and shape type

texts:
I have several shapefiles I need to merge. I have been successful in merging all of them except for one using the code bellow. For the one that is not working, I get an error message that reads:

Traceback (most recent call last):
  File "H:\scripts\sanitary_merge.py", line 19, in <module>
    arcpy.Merge_management(feature_classes, outshp)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\management.py",
 line 4124, in Merge
    raise e
ExecuteError: ERROR 000302: The input feature types do not match
Failed to execute (Merge).

I added in the Walk function, hoping that will fix it, but it is not working still.

Code:

import arcpy
import fnmatch, os

workspace = 'E:\GPS\GPSDATA'
outdir = 'E:\GPS\GPSDATA\GPS_merge_layers'
search = 'Sanitary.shp'
outname = 'Merge_18.shp'

feature_classes = [] #Empty list to collect shapes.

walk = arcpy.da.Walk(workspace, datatype="FeatureClass", type="Point")

for dirpath, dirnames, filenames in walk:
    for filename in filenames:
        feature_classes.append(os.path.join(dirpath, filename))

if feature_classes:
    outshp = os.path.join(outdir, os.path.basename(outname[0][:-4]) + 'tanya_Sanitary.shp')
    arcpy.Merge_management(feature_classes, outshp)


texts: Missing styles in symbol selector - ArcMap 10.3

texts:
I am trying to add "Transportation" style reference to change the display of my point symbols, but for some reasons I can't. When I am in Symbol Selector I select "Style References" and then check in "Transportation", but it doesn't add new symbols to choose from.

When I go to Style Manager the "Transportation.style" icon is grayed out, but at the same time I can go into "Marker Symbols" and I can see that symbols are there.

Please see attached picture.

Am I missing some files? If yes, where could I download them?


texts: Reading 1/3 arcsecond DEM data using C++?

texts:
I am writing my own application to use and present elevation data. So far, I've been easily using the three arc second SRTM data presented in ArcInfo ASCII format. But now I'm looking to use the higher resolution 1/3 arc second data. And I find them downloadable in three formats: arcgrid, gridfloat, and img. No ASCII format.

I looked into the arcgrid format, and it looks utterly opaque and complicated. So here's the question if anyone knows: which of these three would be easiest to code for, to read the elevation data into a C or C++ array in my program. And if anyone has any code snippets to share, I'd love to see them.


texts: Setting junction properties in ArcGIS Desktop?

texts:
I have a situation where I have a line feature class and need to build a network out of it. The lines represent roads.

I wonder how I can model the intersections on the network. Almost all of them are going to be ordinary street intersections; someone coming down one street can either go straight or turn either way at the crossroads. However, a few of the streets are one-way, and there will be a few places where roads cross but do not connect to each other (for example, a street crossing an expressway on a bridge). How do I go about setting this property for an intersection?

I understand that, in ArcGIS 10.2, I need to have a point feature class for junctions (correct me if I'm wrong). Is there a way to create such a class without manually clicking on every single intersection? If there is a way, could I perhaps define an ordinary intersection (i.e., corner of two-way streets) as the default type and then manually define exceptions? That could be more manageable than manually digitizing every place where streets cross.


texts: Supposing Google uses GTFS format for Directions service, how does it represent the actual bus route?

texts:
So I am actually taking a look right now at Google's GTFS format for public means of transport. But from my point of view, I can see information included about the position of the bus stops and some generic info about the "trips" as they call them, but how do they get to represent the right route that connects two stops on Google Maps?

Maybe shapes.txt is used for that? Is it though?


texts: Splitting 3D Spheres with 3D Walls (lines converted to 3D) in ArcGIS for Desktop?

texts:
I am using ArcGIS 10.3 for Desktop.

I would like to split spheres (created from 3D buffers) in half using 3D fences/walls (lines converted to 3D).  Both datasets are multipatches.

The purpose of doing this step is to create hemispherical features that represent downwind directions from air quality monitoring sites. The basic idea is to measure a variety of land use and other variables that may explain variation in observed air quality.

I figured out one workaround by creating 3D boxes (see 2nd fig) from one-sided buffers of lines that bisect the circular buffers (2D version of spheres), creating 3d boxes from these buffers, then using 3d intersect but if there is a more elegant solution - i.e., a way to split multipatch features, that would be helpful.


texts: Issues when using pgrouting's shortest path

texts: I am having some issues when using pgrouting's shortest path. I am following this tutorial available at http://workshop.pgrouting.org/chapters/shortest_path.html#dijkstra. This query returns a list of edges. Then I get the start and end points of each edge. Then I just join these points to show the shortest path. Just to let you know that I don't necessarily get the ending point of the first edge in the list to be the starting point of the second edge in the list returned. So I join, the ending and starting points of the two edges as well. However, sometimes the points are repeated and I have lines crossing over each other and the path looks a mess. What could be the reason behind this?


texts:
In QGIS, the option to save the style to a database (postgres) does not appear in the drop down list if the postgis table has been opened using add vector layer and then entering the postgres database details, rather than add postgis layer

I see that the connection string is different (PG:dbname) but why does this make a difference?

I am using 2.14


texts: Arcpy SearchCursor reads csv field incorrectly

texts:
My code is:

with arcpy.da.SearchCursor("txtFilePath", ["Field1", "lotsOFOthers"]) as cursor:
  for row in cursor:
    val = "%s" % row[0]
    ... non related code which selects a few values
    newTextFile.write(val + "\n")

My issue is that field 1 is a mix of integer and string values, and sometimes the values that are strings get written out as nulls.  There is literally nothing written out for that value.  I strongly suspect that this is due to Arc seeing only integer values in the first few records and surmising that the field is type int, and then when it tries to read a string it fails.

My question is, is there a way to make sure this won't fail?  I know from experimentation that this works when some or all of the first ten records contain strings; in that case both integer values and string values are written out correctly.

I understand I could read and write csv files without using cursors, but is there a way to do this while still using da.SearchCursor?


texts: Can I extract the `highway` layer using R's sf-package?

texts:
From here
it seems this is not possible:

  osmdata returns all ways, and thus enables, for example, examination
  of the full attributes of any  member of a multigeometry object. This
  is not possible with the GDAL/sf translation. As for points,  the
  only additional difference between osmdata and GDAL/sf is that
  osmdata retains all key-value  pairs, whereas GDAL retains only a
  select few.

Which packages can select highways? I found an answer here, but this is not working in my case

In a first step, I want to display the highways on a map...

Edit: It seems I get closer - almost without using library(osmdata)...

library(sf) # Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
sf::st_layers("myfile.osm", options = c("USE_CUSTOM_INDEXING=NO"))

options:        USE_CUSTOM_INDEXING=NO
Driver: OSM
Available layers:
        layer_name       geometry_type features fields
1           points               Point       NA     10
2            lines         Line String       NA      9
3 multilinestrings   Multi Line String       NA      4
4    multipolygons       Multi Polygon       NA     25
5  other_relations Geometry Collection       NA      4

and

ims <- sf::st_read("myfile.osm", options = c("USE_CUSTOM_INDEXING=NO"), layer = "lines")
options:        USE_CUSTOM_INDEXING=NO
Reading layer `lines' from data source `C:\myfile.osm' using driver `OSM'
Simple feature collection with 47 features and 9 fields
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 7.554 ymin: 50.20 xmax: 7.564 ymax: 50.2
geographic CRS: WGS 84

Then

ims <- sf::st_transform(ims, crs = 4326)
my_bbox <- sf::st_bbox(ims)
ims <- ims %>% osmdata::opq(my_bbox) %>% sf::add_osm_features(key = "highway")
bb_poly has more than one polygon; the first will be selected.
 Error in bbox[, 2] : Indixing out of bounds

The solution

ims <- sf::st_read("myfile.osm", options = c("USE_CUSTOM_INDEXING=NO"),
                   layer = "points") # Those are the points
ims <- sf::st_read("myfile.osm", options = c("USE_CUSTOM_INDEXING=NO"),
                   layer = "lines") # That is the highway!!
# ims <- sf::st_read("myfile.osm", options = c("USE_CUSTOM_INDEXING=NO"),
#                    layer = "multilinestrings") # does not work
ims <- sf::st_read("myfile.osm", options = c("USE_CUSTOM_INDEXING=NO"),
                   layer = "multipolygons")
# ims <- sf::st_read("myfile", options = c("USE_CUSTOM_INDEXING=NO"),
#                    layer = "other_relations") # does not work

ims <- sf::st_transform(ims, crs = 4326)
ims_geom <- sf::st_geometry(ims)
plot(st_geometry(ims_geom))


texts: Seeking shapefile of hexagons that covers world?

texts: I'd like a persistent .shp that covers the world in hexagons (preferably of a decently small area)–does anybody know where I could find something like this, or perhaps what software I could use to make one? MMQGIS freezes whenever I try to export something of this size.


texts: Offline map with less zoom size and maximum area covered

texts:
I am developing a mobile app for which an offline map is required. Is it possible to have such map within few 10s of MBs? The requirement is to have:

Map with smaller zoom size; e.g. In this URL of Google Maps, the
z denotes the zoom size (following longitude, latitude): https://www.google.co.in/maps/@12.9677753,77.6490654,15z; Hence 15
Google Map equivalent zoom size is manageable, I am ok with 14 as
well if it reduces the size considerably
The map can be a big image file with names of the street etc.
printed
The map should cover at least a country as big as US/Canada (if not
the whole world)

Google maps doesn't seem to help. I also tried OpenStreetMap and its various links but what I found was the map size of 30-40 GB. I referred following question as well, but unable to get answer:
Offline Map Service with Mapnik and Leaflet

From where can I download such small size map or what is the minimal size of map which has bare minimal street view (not sure if this link can help)?


texts:
I was trying to subtract the two raster images using the code:
import rasterio as rio
import numpy as np

dataset = rio.open("FloodArea_20191122T001936.tif")
data = rio.open("FloodArea_20200929T001950.tif")
subtract = dataset.read(1) - data.read(1)

with rio.open('E:\lidar_chm.tiff', 'w', **dsm_meta) as ff:
     ff.write(subtract,1)

but it shows the ValueError: operands could not be broadcast together with shapes (10963,11058) (10781,11181).


texts:
I do a lot of parcel editing and have it in a topology set to a UTM coordinate system. 99% of all survey plans I add into the parcel fabric are using metres as there unit. That 1% of the time I get a survey plan in feet always takes me 10x longer to input. This is basically just because I cannot on the fly set the editing units from metres to feet.

I really hope that there is an easy way to set the edit units from metres to feet, but I haven't been able to find anything so far. Right now I have to copy all the units from the plan into excel, convert them from feet to metres and then start imputing the plan. I figured it'd be easy, like there would be an on the fly conversion tool within ArcMap, but I haven't found one. I also though it'd be easy to simply change the map units, but they are greyed out:

I can't fathom that I am the only one having this issue and I know ESRI has to have some sort of way of optimizing this process. So please show me where the on the fly editor conversion is.


texts:
I am working on a Python script which will give me the x,y and depth coordinates in QGIS. I have depth as an attribute of the raster layer. My objective is to get a status bar which shows x,y and the depth coordinate when the mouse is at a certain point of the canvas. The code I have used is:

from qgis.gui import *
from PyQt4.QtGui import QAction, QMainWindow, QApplication
from PyQt4.QtCore import SIGNAL, Qt

class MyWnd(QMainWindow):
  def __init__(self,layer):
    QMainWindow.__init__(self)
    self.statusBar().showMessage(u"x:-,y:-")
    self.canvas = QgsMapCanvas()
    self.canvas.setCanvasColor(Qt.white)

    self.canvas.setExtent(layer.extent())
    self.canvas.setLayerSet([QgsMapCanvasLayer(layer)])

    self.setCentralWidget(self.canvas)

    actionZoomIn = QAction("Zoom in", self)
    actionZoomOut = QAction("Zoom out", self)
    actionPan = QAction("Pan", self)

    actionZoomIn.setCheckable(True)
    actionZoomOut.setCheckable(True)
    actionPan.setCheckable(True)

    self.connect(actionZoomIn, SIGNAL("triggered()"), self.zoomIn)
    self.connect(actionZoomOut, SIGNAL("triggered()"), self.zoomOut)
    self.connect(actionPan, SIGNAL("triggered()"), self.pan)

    self.toolbar = self.addToolBar("Canvas actions")
    self.toolbar.addAction(actionZoomIn)
    self.toolbar.addAction(actionZoomOut)
    self.toolbar.addAction(actionPan)

    # create the map tools
    self.toolPan = QgsMapToolPan(self.canvas)
    self.toolPan.setAction(actionPan)
    self.toolZoomIn = QgsMapToolZoom(self.canvas, False) # false = in
    self.toolZoomIn.setAction(actionZoomIn)
    self.toolZoomOut = QgsMapToolZoom(self.canvas, True) # true = out
    self.toolZoomOut.setAction(actionZoomOut)

    self.pan()

  def zoomIn(self):
    self.canvas.setMapTool(self.toolZoomIn)

  def zoomOut(self):
    self.canvas.setMapTool(self.toolZoomOut)

  def pan(self):
    self.canvas.setMapTool(self.toolPan)

  def eventFilter(self,source,event):
    if event.type()==Qt.Event.MouseMove:
      if event.buttons()==Qt.NoButton:
        pos=event.pos()
        x=pos.x()
        y=pos.y()
        p=selfcanvas.getCoordinateTransform().toMapCoordinates(x,y)
        self.statusBar().showMessage(u"x:%s,y:%s" %(p.x(),p.y()))
      else:
        pass
    return QMainWindow.eventFilter(self,source,event)

This file I have saved as mywnd.py in the Python27 folder
Then in Python console in QGIS I use:

import mywnd

w=mywnd.MyWnd(qgis.utils.iface.activeLayer())
w.show()

This gives me my active layer but I still can't figure out how to get the values of x and y in the statusbar. I also am not very good at Python.


texts: Save QGIS project to GeoPackage using PyQGIS

texts:
How it is possible save qgis project to GeoPackage in PyQGIS?

the image is user interface but I need in code


texts: Trouble with ArcGIS Server 10.2 web adaptor when Amazon Image start/stop

texts: Every time we stop our ArcGIS Server 10.2 Amazon Image we need to reconfigure the web adaptor.  Anybody have the same trouble or know of a solution?


texts: How to add reclassified layer to another layer in Google Earth Engine

texts:
Ok so let's just jump right to the script
var indon = ee.FeatureCollection('users/sdwijati/IDN_adm2')
var Sesar = ee.FeatureCollection('users/sdwijati/Sulawesi_sesar')
var nama_prov = ['Sulawesi Tengah', 'Sulawesi Utara',
                 'Sulawesi Tenggara', 'Sulawesi Barat',
                 'Gorontalo', 'Sulawesi Selatan']
var prov = indon.filter(ee.Filter.inList('NAME_1', nama_prov));
var s = ee.Image("OpenLandMap/SOL/SOL_SAND-WFRACTION_USDA-3A1A1A_M/v02")
var savg = s.reduce(ee.Reducer.mean());                                   //Mean soil sand content
var smap = savg.divide(20).ceil();
var resmap = smap.remap([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20],
                        [0,0,0,0,0,0,0,0,0,45,55,60,65,70,75,80,85,90,95,100]); //Reclassified mean sand content
Map.addLayer(savg.clip(prov), null, 'Soil sand content', false)

var gtsavg = savg.updateMask(savg.gte(45))                         //sand content greater than 45%
var sandvis= {min: 0, max: 100, palette: ['fbff00','000000']};

var n = ee.Image("OpenLandMap/SOL/SOL_WATERCONTENT-33KPA_USDA-4B1C_M/v01")
var navg = n.reduce(ee.Reducer.mean())                            //soil water content
Map.addLayer(navg.clip(prov), null, 'Soil water content', false)

var xx = resmap.add(navg) //Soil water content percentage + Reclassified soil sand content percentage
Map.addLayer(xx.clip(prov), null, 'Sand + Water Content', false)

var sesar = Sesar
var bufferBy = function(size) {
  return function(feature) {
    return feature.buffer(size);
  };
};

var intervals =
  '<RasterSymbolizer>' +
    '<ColorMap type="intervals" extended="false" >' +
      '<ColorMapEntry color="#0CE84A" quantity="41.9925" label="41.9925" />' +
      '<ColorMapEntry color="#FFFF40" quantity="78.9925 " label="78.9925 " />' +
      '<ColorMapEntry color="#FFBF40" quantity="115.985" label="115.985" />' +
      '<ColorMapEntry color="#EC4040" quantity="152.97" label="152.97" />' +
    '</ColorMap>' +
  '</RasterSymbolizer>';

var sarbuff = sesar.map(bufferBy(25000*4))
Map.addLayer(xx.sldStyle(intervals).clip(prov), {}, 'Susceptibility Categorization');
Map.addLayer(gtsavg.clip(prov), sandvis, 'Soil with sand content >45%')
Map.centerObject(prov, 7)

So, if you run the script it will show you two layers of active map. 'Soil with sand content >45%' and 'Susceptibility Categorization'.
Based on var xx = resmap.add(navg), when I clicked the Google Earth Engine's inspector in the darker area (area with soil sand content >45%), I expected the Soil sand content percentage to be added to Soil water content percentage, which the total percentage will be shown in Sand + Water Content, and thus the total percentage then will be matched to the right color, accordingly, based on the colors on these lines:
'<ColorMapEntry color="#0CE84A" quantity="41.9925" label="41.9925" />' +
'<ColorMapEntry color="#FFFF40" quantity="78.9925 " label="78.9925 " />' +
'<ColorMapEntry color="#FFBF40" quantity="115.985" label="115.985" />' +
'<ColorMapEntry color="#EC4040" quantity="152.97" label="152.97" />' +

Unfortunately, that didn't happen. The Sand + Water Content will be the same percentage with Soil water content, despite the Soil sand content percentage is above 45%. So, what did I do wrong? And, any idea to fix it?


texts: Does Leaflet actually ask for images in tile shapes when adding WMS layers

texts:
I am learning to use Leaflet. I tried adding a WMS layer of schools from GeoServer like this:

 L.tileLayer.wms("http://localhost:8080/geoserver/geog585/wms?", {
        layers: 'geog585:TaipeiSchoolUTF8',
        format: 'image/png',
        transparent: true,
    }).addTo(map);

When I look at the network tab in the developer's tool, I saw the the number of request sent and the number of png files returned from GeoServer is the same as the number of png files returned from the tile base map source I'm using, and the size of these school images returned are what I would expect from OSM/Google Map tiles (256*256). This suggests that the school layer may be returning in the form of tiles. But from the developer's tool, I was also able to see that the requests sent are in the from of WMS request. Here's a sample

Request URL:http://localhost:8080/geoserver/geog585/wms?&service=WMS&request=GetMap&layers=geog585%3ATaipeiSchoolUTF8&styles=&format=image%2Fpng&transparent=true&version=1.1.1&height=256&width=256&srs=EPSG%3A3857&bbox=13521404.55553454,2876478.248427754,13531188.495155042,2886262.188048257
Request Method:GET

So now I am a bit confused. In my understanding, WMS server takes the bbox, crs/srs, height, width requested (and layer and style), draws an image based on these, and sends it to the client. So one might think that Leaflet would sent a WMS request with the bbox of the map view, and width, height, and crs/srs of the map parameters and get a single image in return.

Is it the case that when adding layers to Leaflet, Leaflet want to put everything into tile forms? And so that L.tileLayer.wms() actually sends multiple request that matches the bbox of the tiles the map area covers, and using the conventional width and height of OSM/Google Maps tiles?


texts: gdal_calc.py raster reclass equivalent of ArcGIS reclass or r.reclass in grass

texts:
Would like to achieve this in gdal_calc.Py if possible.

I have a raster with a value of 0-7 in GTiff

I would like one output raster where

NoData = 0,1,2,3

1 = 4,5,6,7

So basically an ArcGIS reclass but using gdal_calc.py. I know how to do it in QGIS, ArcGIS, etc.

I know this question has been asked before for SAGA, QGIS, GRASS (r.reclass), but I do not see it for gdal_calc.py


texts: Getting started with SharpMap in Visual Studio with C#

texts:
I am trying to use SharpMap in Visual Studio with C#.
My first step was to go to the Tutorials on CodePlex.
https://archive.codeplex.com/wikipage?title=SharpMap%20Tutorial&referringTitle=Home
But that is a broken link, as are all the links on the CodePlex page for SharpMap.
I found another tutorial.  The first step to make a map window worked well, but the second step:
https://github.com/SharpMap/SharpMap/wiki/Adding-a-layer-to-the-Map
failed to produce any image in the map window.  (I carefully followed the instructions and do have the data file extracted and correctly referenced.)
I have no idea what to try next.
Is there any documentation for SharpMap?  I've not been able to find any.
My application is to display of the NOAA Nautical charts.


texts: Intersecting two shapefiles from different folders based on their name using scripting in ArcGIS Desktop?

texts:
How can i script for intersecting the two shape files from two different folders based on their names?

i have two folders which contains the shape files with same name i want to intersect the shapefiles based on their name, scripting should pick the files based on their name from two different folders and output file should store in another folder.

I'm using ArcGIS Desktop 9.3 and arcgisscripting (that preceded ArcPy).


texts: Loading an unknown text vector format with QGIS?

texts:
I have a vector data in text format like below, and i want to use it under QGIS as a vector layer.

Each polygon starts with one line and contains: number, type, number of points. Then follows the xys. It is not difficult to parse it with python.

Another text file contains the height of each building.

This file can't be imported as delimited text file into QGIS, and I can't figure out how to create a self-defined vector format.

Is there any way to load geometries from this kind of file into QGIS ?


texts: Uploading Polygon Shapefile to Garmin GPS. How to View all Polygons?

texts:
I am trying to upload a polygon shapefile to a Garmin 60CSx unit. The polygon shapefile is a small area containing polygons of vegetation types. There are are 53 polygons (records) in the shapefile. I am using DNRGPS v6.0.0.15

Only one of the 53 my polygons I have uploaded appears in the GPS unit. I do not know why

This is what i have done

Open DNRGPS and connect to Garmin.
Go to File-->Load From-->File... Locate Polygon Shapefile.
When I choose the shapefile, the 'Select TIDENT Field' dialogue prompts and I choose a field called 'SITE'. The SITE field has the same value for all 53 records ('Forest Stand 101').
Now the polygons are in a table under 'Tracks'. There's 2150 Tracks. I choose Tracks-->Upload. The 'Upload To' dialogue prompts and I choose 'Named tracks'
GPS beeps and I have a successful upload.
I open GPS, go to the 'Tracks' menu and see my new track 'Forest Stand 101'. I open that track and choose the 'Map' button to view my polygons in the map.
This is where I notice only one of the 53 polygons is displaying. How to I get all polygons to display?


texts: Remove features inside a feature in shapefile using python (UNION, CONTAINS methods)

texts:
I have created a shapefile from raster. The shapefile was generated using contours detection. I am in the process of automating feature detection using Python and open source libraries such as shapely, fiona, ogr. The problem I am having is that the shapefile generated works well in detecting features small in size but when it tries to detect large features it also detects features inside the large feature. This results in feature inside a feature problem.

I have been trying to implement polygon inside a polygon implementation but have not been able to implement it successfully. I am thinking of finding a way in which this issue of feature inside a feature gets removed during the process of making the shapefile at the very beginning (raster to shapefile conversion). Need some suggestions on how to proceed on this issue.

The red highlighted area is the region I am interested in and everything inside it needs to be removed.


texts: Is it possible to georeference in CartoDB using just UK postal codes?

texts: Is it possible to georeference using just UK postal codes. I have the minimum free account and have had now luck with this. IF so should the postal codes be continuous strings with no gaps and is it correct to add in the free text box the country as "United Kingdom"?


texts: I have two vector (shapefile) data. One contains soil order for county, and the other one contains boundary for each county. I want to extract the soil order for each county by using county boundary. But I am not sure how to do this?


texts: How to maintain GlobalIDs while copying shapefile data into ESRI FGDB feature class?

texts:
I want to copy a shapefile data into an ESRI file geodatabase feature class. That feature class having a GlobalID column and it is being auto generated by the system while copying.

The GlobalID values maintained as string in the shapefile and I want to retain the same GlobalIDs in copied feature class too.


texts: Converting from different projections to WGS84 Web Mercator Auxillary Sphere

texts: I have multiple datasets in my ArcGIS, each having different map projections. All of them are vectors. Some are using Sphere Mercator, others using WGS84. I cannot perform the reprojection function in ArcGIS as I do not have the formula for geographic transformation. I am stuck and cannot proceed to make a map due to this issue.


texts: Creating Self Executing Python Script (*.py)?

texts:
Trying to figure out how to create self executing python scripts with using the ArcGIS python window.     I'm ok  (generally) when developing and running scripts inside the python ArcGIS windows but unable to figure out self executing scripts.

For example ....here is a really simple script created to run a "Select Case" to change annotate attributes in one field (BBB) based upon attributes in field AAA.  Goal is simply to double click the *.py file and have it run - outside of ArcGIS

# Import arcpy module
import arcpy

# Import env from arcpy as ENV and set the workspace environment
from arcpy import env
env.workspace = "C:/Users/Jeff Reichman/Documents/My GIS Projects/OpenStreetMap/OSM.gdb"

# Local variables:
Input_Field = "C:/Users/Jeff Reichman/Documents/My GIS Projects/OpenStreetMap/OSM.gdb/SaintLouis/SaintLouis_osm_ln"

# Process: Calculate Field
arcpy.CalculateField_management(Input_Field, "BBB", "x", "VB", "dim x as long/nselect case [AAA]/n  case 1: x = 100/n  case 2: x = 200/n  case 3: x = 300/n  case 4: x = 400/n  case 5: x = 500/n  case else: x = -999999/nend select")


texts: Display world map raster in World Robinson projection in QGIS without the corners being repeated?

texts:
I am using QGIS version 3.2.3-Bonn on Windows 10.

I have downloaded the global physical map raster from Natural Earth link. This is in WGS84, and I wish to project it as World Robinson (ESRI:54030), but when I do so the map repeats in the corners (see screenshot). I wish to create a PDF without the corners being repeated. How can I do this?

I have tried:
 - Changing the project CRS to World Robinson and enabling 'On The Fly' projection.
 - Reprojecting the file to the World Robinson CRS, using raster -> projections -> warp (reproject).

Both approaches result in the corners of the map being repeated.

How can I simply show the map reprojected, without any part of the map being repeated?


texts: How do you find a lat/lon (DMS) on a paper map using a 30cm ruler?

texts:
How do you use a 30cm ruler to find DMS on a paper map? The locations I would like to find are the 'corner' points so I can generate an extent based on the four corners.

I have an old paper map (3 actually) for Northern Canada (late 1800's) that doesn't provide the Ellipsoid or the Datum. It provides a representative fraction (1:660,000 approx) and a scalebar (1" = 10 2/3 mile). The map shows grid lines spaced every 1 degree. No minutes or seconds are labelled.

I understand that NOT knowing the datum or ellipsoid will automatically introduce a margin of error into the calculations, but this is not a big deal for this exercise.

I determined the Lat/Lon of the intersecting grid lines, and from this question, was able to infer that it is closest to Lambert Conformal Conic (Statistics Canada, EPSG 3347).

Below is the index map showing all 3 maps with grid lines every 2 degrees:

I will need to do this process for all three maps as those grid lines are spaced every 1 degree, and not 2 as in the index above.

Of course, I could geo-reference to a known spatial reference in a computerized GIS and digitize the extent, but what if your GIS is PC-less and you have traveled back in time and are now stuck...

If it is easier to provide an answer using say, an engineers ruler (1:100, 1:2500 etc) then feel free. It's just a 30cm ruler seems to be more readily available in a given situation.


texts: QGIS 3.4 : Add points to exact locations?

texts: I would like to add several coordinates at specific British National Grid Easting and Northing values. I know this can be done by bringing in a .csv, but I want to know is there a way to manually enter the numbers in QGIS 3.4, and then the point moves to that location ? I've read that in QGIS 2.x there was a plugin called Numerical Vertex Editor that doesn't seem to exist for 3.x ?


texts: What options are available for the rasterio .update() method?

texts:
I have a script which clips raster data by a shapefile using rasterio and fiona:

import fiona
import rasterio
import rasterio.mask

def clip_raster(shp, imagery, out_imagery):
    with fiona.open(shp, "r") as shapefile:
        features = [feature["geometry"] for feature in shapefile]

    with rasterio.open(imagery) as src:
        out_image, out_transform = rasterio.mask.mask(src, features, crop=True)
        out_meta = src.meta.copy()

    out_meta.update({"driver": "GTiff",
                     "height": out_image.shape[1],
                     "width": out_image.shape[2],
                     "transform": out_transform,
                     "compress": "LZW",
                     "nodata": 255})
    with rasterio.open(out_imagery, "w", **out_meta) as dest:
        dest.write(out_image)

if __name__ == "__main__":
    shp = r'/path/to/shapefile.shp'
    imagery = r'/path/to/input_imagery.tif'
    out_imagery = r'/path/to/output_imagery.tif'

    clip_raster(shp, imagery, out_imagery)

You can see that the .update() method has a variety of options such as setting the compression type ("compress": "LZW") or nodata value ("nodata": 255):

out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform,
                 "compress": "LZW",
                 "nodata": 255})

Is there a list or documentation of all the available options available to the .update() method when writing raster data using rasterio?


texts: How to download global SOC and clay soil datasets?

texts: I just wanted to clarify the best way to access the latest release of SoilGrids (5 Dec 2019). Are the datasets for clay (CLYPPT)at: https://files.isric.org/soilgrids/data/recent/ the same as https://files.isric.org/soilgrids/data/preview/clay/? Or are these the older versions?


texts: QgsComposerItem error when runnning on Debian instead of windows

texts:
I have the following script which has worked on Ubuntu 16.04 and Windows but not on Debian.

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtSql import *
from PyQt4.QtXml import QDomDocument
from qgis.core import *
from qgis.gui import *

qpt_file = "a4_vertical.qpt"
output_html = "test"
if QGis.QGIS_VERSION_INT > 20400:
    renderer = iface.mapCanvas().mapSettings()
else:
    renderer = iface.mapCanvas().mapRenderer()

c = QgsComposition(renderer)

qpt_file_dom = QDomDocument()
qpt_file_dom.setContent(QFile(qpt_file), False)
c.loadFromTemplate(qpt_file_dom, None)
c.refreshItems()

composerMap = c.getComposerMapById(0)

composerLabel = c.getComposerItemById('HEADER')
composerLabel.setText(output_html)

When running the code on Debian I get the following error:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/tmp/tmpB8icUm.py", line 26, in <module>
    composerLabel.setText(output_html)
AttributeError: 'QgsComposerItem' object has no attribute 'setText'

I tried to solve this based on this issue which is 6 years old:
https://issues.qgis.org/issues/3236

It works when I add composerLabel.__class__ = QgsComposerLabel.

But when I'm trying this with QgsComposerFrame as follows:

table_frame = comp.getComposerItemById('HTML_FRAME_1')
if not isinstance(table_frame, QgsComposerFrame):
    table_frame.__class__ = QgsComposerFrame
table_frame = table_frame.multiFrame()

I get the AttributeError:

AttributeError: 'QgsComposerFrame' object has no attribute 'multiFrame'

I'm running Qgis 2.14.18-Essen on Debian Stretch.

QGIS version information:

QGIS version: 2.14.18-Essen
QGIS code branch: Release 2.14
Compiled against Qt: 4.8.7
Running against Qt: 4.8.7
Compiled against GDAL/OGR: 2.1.2
Running against GDAL/OGR: 2.1.2
Compiled against GEOS: 3.5.1-CAPI-1.9.1
Running against GEOS: 3.5.1-CAPI-1.9.1 r4246
PostgreSQL Client Version: 9.6.3
SpatiaLite Version: 4.3.0a
QWT Version: 6.1.2
PROJ.4 Version: 493
QScintilla2 Version: 2.9.3

I tried run this code under Ubuntu 16.04 and Windows 7 and there it runs without problem.


texts: Slightly confused. I need to add 3 point features to an already created map project (state plane projection). The description given to me about the location of the points is the planar location, so no matter what I try and use as the display unit, nothing will work for finding planar x,y coords. Do I need to create a new map project with an Azmiuthal-Equadistant projection and create the points there, then import them into the premade map or am I missing something very simple?


texts: WFS point layer does not appear on OpenLayers map when called from GeoServer

texts:
I'm writing because I attempted almost all the tries suggested in any other Q/A on this thread. I have this vector points layer on my local GeoServer and i'm trying to plot my points on a OSM base map through OpenLayers. My script (very basic) is the following:

<!DOCTYPE html>

    <title>Wfs Map</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v5.1.3/css/ol.css" type="text/css">

Everything that I obtain is the OSM layer and no points.
The original sdr of the layer is 4326.

The strange thing (unexplicble to me) is that if I place the same .json file in my github sandbox and I use:

url: 'https://raw.githubusercontent.com/annalisapg/sandbox/master/barriere2_wfs.json',

the points appear. So, in your opinion, what's wrong between the script, the layer and my geoserver installation?


texts:
I am a new QGIS 2.0 user. I was so impressed the new default tool, Topology Checker, in this new version.

But, hundreds errors are bothering me for cleaning up features. Is it possible to clean errors automatically, not manually?

I have an additional 2 questions.

Is it possible to clean topological error in QGIS not GRASS?
According to "v.clean tool", is it possible to clean overlapped areas? In my understanding, there is no specific tool for overlap error.


texts: How to store Affine Transfroms with PostGis 2.3

texts:
I need to store a specific 3D affine transformation for each entry of a table.
Such transformations are later used in order to provide representation of geometries in the entry local reference frame. Notice that I might have a lot of these transformations and thus I would rather avoid defining a custom SRID for each one.

I know that I can apply a transformation using the function ST_Affine
and providing all the transformation parameters.

What is the best way of storing the 6 affine parameters? I do not want to use six columns. I would rather prefer having a single column containing all the transformation parameters. Is there a default type to represent this?


texts:
I'm working on project that uses Open Layers. I have a problem with styling some features. Is it possible to set as Style image and also specify fill and stroke for it, for example background of the png image? At the moment I'm trying to do this like this:

this.styleHomeGreen = new Style({
      image: new Icon({
        src: 'content/map-icons/home-green.png',
      }),
      fill: new Fill({
        color: 'red'
      })
    });

In result image is being shown properly but there's no fill. I also tried with Stroke but with same result.


texts: SQL Server 2014 Failure to access the DBMS server 08S01

texts:
After upgrading our esri sde database from SQL Server 2008 to 2014 I am experiencing issues when overwriting some of our larger datasets. After some amount of time the overwriting process abruptly terminates. The error seems to be coming from the SQL side and reads: "Failure to access the DBMS server [08S01:[Microsoft][SQL Server Native Client 11.0]TCP Provider: The specified network name is no longer available.]"

This happens both when using python or the Arc tools user interface to perform either the Append or FeatureClassToFeatureClass tools. I've also tried running these tasks on off hours when the server is under minimal load to no avail. There were no issues when we were on SQL Server 2008 and the python scripts are the exact same as before the upgrade.

I have the SQL Server Native Client 2012 installed and the Remote Query Timeout setting is at 600 seconds. Our ArcGIS client is 10.3

Has anyone experienced and resolved this issue before?


texts: Isolating just one area of larger map in QGIS?

texts:
I have a map of Wales with various layers on it in QGIS 2.18.20. I want to select just a small rectangle of the country, with the layer details which apply to that area, and then export it.

I've tried the multiple selection tool but it doesn't seem to be possible to save the selected area.


texts: Selecting first 5 rows of selection after reverse sorting using attribute field in ArcPy?

texts:
I have a feature class containing 380 polygons representing sedimentary basins and another polygon feature class representing oil discoveries within said basins. My plan is to:

Iterate through all basins and select all fields that fall within them.
Reverse sort by a oil recoverable volume attribute (largest to smallest).
Sub-select the first 5 rows of the selection and export this to a feature class (other rows will be appended).

The first two steps I think I am ok with however, I'm stuck on how I would go about creating the sub-set where the first 5 rows are selected and then exported. If there are <5 rows, all will be selected, if there are 0 rows they will be ignored.

How, using ArcGIS Desktop 10.2, can I achieve the goal of selecting the first 5 rows of a selection?

The da.SearchCursor can be used to sort but not sure after that.


texts: Interactive geographic region definitions

texts:
Background

An organization has regions within a province, but has no list of coordinates that define region boundaries. I would like to draw custom polygon areas on a map using a web-based and user-friendly interface (such as Google My Maps or Scribble Maps). The regions are for data aggregation.

Technical

Maptitude and OpenGeo are being considered for the final GIS application. The database will likely be PostgreSQL running PostGIS.

Question

What open-source tool would you use to define the coordinates by drawing polygons on an interactive map? That is:

Browser -> Map -> Polygons (shared boundaries) -> Export

If the export step (to PostGIS) was performed automatically, so much the better.

A cross-browser solution (IE8, FF, Chrome, Safari) with a working demo would be most helpful.

Related Links

https://stackoverflow.com/questions/2325670/draw-polygon-using-mouse-on-google-maps


texts:
I have a tile layer created on ArcGIS for Server 10.3. most of the tiles have worked out fine, but I have a number of tiles that are odd colours. I have tried rebuilding the cache, but the odd tiles are still there.

The aerial in the tiles is four separate mrSid files.

What is causing the odd colours and how do I replace them?


texts: Data Transfer of point on line using the anchor point

texts:
Can I a line feature, the anchor point touches a point feature transferred the field value, even when the point touches another line(s)?
How does this work?


texts: Why do clipped polygon features (clipped with adjacent polygon tiles) still overlap with each other?

texts:
I have two layers in QGIS with intersecting polygons (blue "P1" and green "P2"). I want to clip them with adjacent polygons tiles (dotted line - blue "cut1" and green "cut2").

Why the result clipped polygons ("Clipped1", "Clipped2") still overlap? The red line is the result of the topology checker tool in QGIS with the rule "Clipped1 must not overlap Clipped2".

The checker is indeed right because when I try to collect geometries in "Clipped1" and "Clipped2" to one layer and one feature, the result is invalid self-overlapping multipolygon geometry.
This is just an example. I am doing this because I have vectorized land use data from overlapping raster tiles. I want to clip these data to exact adjacent boundaries (specified in adjacent polygons tiles like "cut1"). I do not want the result data to overlap.
The example QGIS project and data are shared here:


texts: Adding collapsible legend box to Mapbox and features no longer click-able

texts:
So, first I started with this map: https://codepen.io/melissad34/pen/mNGjKR

It works exactly the way I want it to work.

Then I added a collapsible sidebar: https://codepen.io/melissad34/pen/LYYrLKN

And now the trails features aren't clickable as they were on the first version.
I had to change the z-index of the divs to get the sidebar to show over the map div, so I wonder if that's why, but I assume that shouldn't be an issue because it's not like the sidebar covers the entire map.

What is causing this?


texts: Rotating hatch labels based on orientation changes in Data Driven Pages Strip maps?

texts:
I've created a strip map for a linear feature and have it loaded as a data driven page.  I want to show stationing for the line.  The stationing displays properly on the first data driven page but on each subsequent page the stationing labels get messed up as the map orientation changes. I want the stationing labels perpendicular to the linear feature as shown in the first image.  The tick marks seem to rotate and stay perpendicular to the line but the labels won't.

How can I fix this?


texts: using custom font with more than 255 symbols

texts:
I would like to use a font for the styling of a layer.

The font contains more than 255 characters. QGIS seems to read only the first 255 characters of the font, right? At least, I can only see the character from 1 to 255 in the styling-dialog.

How can I use one of the other characters of the font?


texts: Labeling numbers on points in a row close together

texts:
In my map I have 76 points aligned on a row. Now I would like to force label each point with an attribute number. I don't mind how far the label is from the point as long that it is shown but of course the closer the better.

In QGIS I found the option of Manhattan-lines. But for what I saw is that I've to manually tell how far the label is from the point. I don't want that. Any ideas or tips?

I'm looking for something like a dynamic placement of the labels. As close as possible to the points but not overlapping each other.

I have both ArcGIS Pro and QGIS. So any solution in one or the other would be perfect.


texts: Interpolating raster corridors in ArcGIS Desktop?

texts:
I have a set of 24 raster (2.5 x 2.5 km, tif) that displays vegetation (from las). However, there is a 10 m symmetric grid of unprocessed lines between the tiles.

Any way to interpolate these?

I have looked into ArcMap Focal Statistics tool, but I don't think this is a solution. I do want to "copy-paste" 5 m och each edge, move 5 m and finally mosaic these.

Ideas?

Rectify?

The elevation isolines are in shp and a later problem - snap?


texts: Syntax error in a query with ST_AsText and ST_CENTROID


texts:
I am trying to find overlapping lines using ST_EQUALS on two tables

This query returns 0 rows:

SELECT * FROM geom1 c, geom2 h WHERE st_equals(c.simple_geom, h.simple_geom);

This query runs correctly and returns overlapping lines:

SELECT * FROM geom1 c, geom2 h WHERE st_astext(c.simple_geom)=st_astext(h.simple_geom);

I am not sure why ST_equals fails.
Both colums are indexed,
I thought it was due to geom being multilinestring, but after converting to simple_geom with just LineString (using St_geometryN(geom,1)) it is still giving unexpected results.

So it is not indexes, and neither MultiLineString column.

What are other possible sources of this behaviour?


texts: How to batch create filenames in ModelBuilder?

texts:
I am running a model in ModelBuilder. The input are raster files in *.tiff format and the output are *.asc format.

I am using iterate raster and there will be 365 input and output files using %name%.asc

e.g input file : LST01011990.tiff  --> output file LST01011990.tiff.asc

How I can output the filename to LST01011990.asc (without repeating the *.tiff)?


texts: Localization of a point on a line in QGIS

texts: I'm relatively new to the GIS field. I have several polylines (streets) and would like to locate stations on these lines. The start and end points of the lines, the length of the line and the distance of the stations from the starting point are given. I'm looking for the coordinates for the stations. With GRASS and PostGIS there are probably possibilities. Does anyone have good explanations?


texts:
I dissolved a polyline feature (streets) on the following fields:
WHOLESTNAM, NUMBEROFLANES, THOROUGHFA, SPEEDLIMIT, GEOID

and asked the dissolve to take an average of the AADT field.

Why does the result contain features that are exactly the same on the 5 fields I dissolved on?

Shouldn't these have dissolved into a single feature?


texts: I have some address data that's given as UDPRN codes, and want to plot them as points in QGIS. Is there a way to convert these to lat/long or join to some existing dataset that can get me there? The only other address data in the set is postcode, but I want to get more specific than that.


texts: Is it possible to convert map tick into feature class easily?

texts:
I have to design a lot of map into the following layout for forest carbon inventory survey:

--There will be a grid set in the map 30 second by 30 second (these are straight lines crossed)
--There  will exist points (on the grid intersection) symbolized based on the layer land use type which is another polygon layer

Please see image--

    I tried fishnet in both arc map and X tool, but there is no way to set 30 seconds interval exactly-- i think this is error prone

    I tried Exporting and using Arc Scan to get the point-but problem is Arc scan limited to poly-line not point what i need.

    After georeferencing, i tried to use spatial join this point layer with land use layer and symbolized accordingly

Is there any easy way to do this..

Thanks in advance


texts:
In this attribute table, "DISTANCE" field name created with this python expression (in Field Claculator):

str( !Shape_Leng! ) + !field_name!

Everything works fine.

But I seek a way to get on 1 number of decimal place and not 9 places as in the result. In the Field Properties the option to change is blocked.


texts:
When executing the following command string in the active browser, Qgis projects the triangle polygon and associated attribute table correctly.

Upon reopening the shapefile in a new browser, the attribute table is lost and no polygon exists.

Is the map registry update command incorrect prior to the 'writer' command, or is the writer command itself improper?

Thanks for any clues.

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *

# create layer
layer =  QgsVectorLayer('Polygon?crs=epsq:4326&field=BOTTOM:integer(8)&field=TOP:integer(8)&field=K:integer(10)&field=R:integer(10)&index=yes', 'poly_temp' , "memory")
pr = layer.dataProvider()

# start editing
layer.startEditing()

# add the polygon
base = QgsFeature()
fields = layer.dataProvider().fields()
base.setFields(fields)
base["BOTTOM"] = "50"
base["TOP"] = "100"
base["K"] = "0.003"
base["R"] = "0.00005"

# Draw feature
poly1 = [QgsPoint(50,50), QgsPoint(25,30), QgsPoint(60,30)]
base.setGeometry(QgsGeometry.fromPolygon([poly1]))
ok, of = layer.dataProvider().addFeatures([base])
iface.mapCanvas().refresh()

# update fields
layer.updateFeature(base)
layer.updateFields()

# update extent of the layer
layer.updateExtents()

# finish editing
layer.commitChanges()

# add the layer to the canvas
QgsMapLayerRegistry.instance().addMapLayers([layer])

# write layer to file
writer = QgsVectorFileWriter("C:\Users\...\Base_Layer.shp", "CP1250", pr.fields(), QGis.WKBPolygon, None, "ESRI Shapefile" )


texts: Can I add entries to the layerlist context menu?

texts:
I am implementing a python based plugin for use with QGIS 1.8.

Most of it's user interface is contained within a QDockWidget but there are some layer specific controls I would like to append to the context menu displayed when the user right-clicks on the Table Of Contents.

Is there a way for me to manipulate this menu?


texts: Converting LIDAR based measurement to UTM coordinate?

texts:
I have lidar based measurement with azimuth elevation and range, I would like to convert this data to UTM coordinate. I know the position of the lidar in UTM coordinate system, so I developed the following code in Matlab

LOS = [cosd(ele)*sind(azi); cosd(ele)*cosd(azi); sind(ele)]
%

R1 = [cosd(yaw), -sind(yaw), 0 ;...
      sind(yaw), cosd(yaw) , 0;...
      0     , 0      , 1];

R2 = [1,  0     ,  0     ;...
      0,  cosd(pitch),  sind(pitch); ...
      0, -sind(pitch),  cosd(pitch) ];

R3 = [cosd(roll)  , 0 , -sind(roll);...
      0       , 1 , 0     ;...
      sind(roll) , 0 , cosd(roll) ];

LOS2 = [LOS'*RR]

x=LOS2(1)*range+OFF2(1)
y=LOS2(2)*range+OFF2(2)
z=LOS2(3)*range+OFF2(3)

But I am not quiet sure that the rotation matrix is good. Am I gooing in the right path?


texts: I am having a problem using QGIS ver 2.10.1 Pisa.  I am trying to use the polygonize function on a raster .tif image, and it runs, but only produces a thin sliver of polygons at the top, then stops working or loading.  Has anyone encountered this problem, or know how to solve it?


texts:
I want to be able to view data layers from LandFire.  The endgoal will be to be able to query these layers for data as a function of latlng. They offer a myriad of data layers through the arcgis servers, which can be found here.

I am trying to add them into a Leaflet map using esri-leaflet's featurelayer.  My code looks like this:

var LFGroundCoverLayer = L.esri.featureLayer({
  url: 'https://landfire.cr.usgs.gov/arcgis/rest/services/Landfire/US_200/MapServer/15',
  opacity: 0.75,
  useCors: false,
})

LFGroundCoverLayer.addTo(leafletMap)

But I get nothing back.  The Network tab of my browser shows I'm getting a 200 response from their server, but nothing shows up in the map.  (I do see that the network request has a ?f=json at the end of it, which doesn't make sense as these are not JSON layers, they're raster layers.) When I try using an imageMapLayer instead, I get a 400 bad request error.  I'm not sure what I'm doing wrong here - I thought I was following the esri-leaflet instructions for featureLayer - should I be using something else?  I am able to successfully load in other featurelayers or imagemaplayers from ArcGIS Server.

LandFire's data page has other options for requesting data layers - is there another approach I can / should try besides esri-leaflet?

Working codesandbox demonstrating the problem

There are some other featurelayers and imageMapLayers commented out there, so feel free to uncomment them and you'll see that they're working correctly.

Edit - Getting Auth

Thanks to TomazicM's comments, I realized this is an auth issue.  I was able to authenticate an Esri layer using these instructions.  But when I try this same tactic with the landfire server, I get bad responses.  When I try  https://landfire.cr.usgs.gov/arcgis/tokens?client_id=<my-client-id>&client_secret=<my-client-secret>&grant_type=client_credentials, I get a reponse that says HTTP GET is disabled

When I try https://landfire.cr.usgs.gov/arcgis/tokens/generateToken?client_id=<my-client-id>&client_secret=<my-client-secret>&grant_type=client_credentials , I get a reponse that says:

Invalid request
Usage: https://landfire.cr.usgs.gov/arcgis/tokens?request=gettoken&username=username&password=password&
Usage: https://landfire.cr.usgs.gov/arcgis/tokens/generateToken?username=username&password=password&
Usage: https://landfire.cr.usgs.gov/arcgis/tokens/gettoken.html

When I try an API call with that suggested format https://landfire.cr.usgs.gov/arcgis/tokens?request=gettoken&username=username&password=password I get a response that says:

Generate Token request should not contain username and password in the query string.

I'm sort of at a loss here.  Is anyone familiar with using ArcGIS rest API to authenticate users?


texts: Does ArcGIS for Desktop have auto-backup feature?

texts:
Did Esri decide to dump the backup file when ArcMap crashes and didn't save the file? Did they change anything to that? Did they change that to MXD Doctor?

Whenever I try to exit or save my project in ArcMap, it just crashes and I when I look for the backup file, I see nothing.

I'm using ArcGIS 10.1 for Desktop.


texts: qfield will not highlight selected feature

texts: When selecting a feature in qfield it will not highlight the way it did in past versions. I do not know why this is happening. Was this taken out from the past update?


texts: Using Field Calculator of ArcMap to populate field based on codes from another field

texts:
I want to do this basically in the Field Calculator ArcMap using its Python Parser:
if !AELEG! == 'L04'
 !rock_type_age!.replace('Mid-Pleistocene-Holocene volcanic rocks (mainly rhyolite and dacite pyroclastic flow deposits')

So the field name with the codes is AELEG and the code in question here is L04, for every row in the field rock_type_age that corresponds to L04 in the field AELEG I want to add a string (in this case 'Mid-Pleistocene-Holocene volcanic rocks (mainly rhyolite and dacite pyroclastic flow deposits'))
Any suggestions?


texts: Insert Geography data into SQL Server 2008

texts:
Can someone guide me to or give me an example on how would one go about insert/query Geography data in Sql Server 2008 ( Polygon with 4 points  and Point in specific) using C# EF or Linq2Sql.

I've been searching for few days now and i can't figure it out :(

Thanks in advance!


texts: Install ArcGIS apps on my iPhone. Is there a manual or instructions for it?

texts: I have installed ArcGIS ESRI on my iPhone. Is there a manual or steps on how to use it? I want to know how to use it.


texts: No Output When Using Python Console In QGIS

texts:
I am fairly new to using the python console in QGIS so I am not entirely sure if there is an issue with a line of code I am using or human error.

I am using this code to cycle through a  point layer of 3000 points situated over around 14 vector layers that I have created. The aim of the code is to efficiently show what layers each point resides on. As each vector layer is assigned a rank I cannot do a select by location as that would take months.

When running a test of around 10 points, the code runs with no issues, but when I go to open the text file I have created for the output, it is empty.

However, interestingly when I ran the code with both the test data points and with the original dataset of 3000 still active in QGIS, an output was produced. So I am assuming the code works and produces an output when a point overlays another point but not when a point overlays a vector layer. Is this an mistake within the code or might it be another factor? For example, the text file has to be formatted in a certain way? The code is displayed below.

counter = 0

# Set the path for the output file
output_file = open(r'/Users/admin/Desktop/TestFlood.txt', 'wb')
iter_layer = QgsProject.instance().mapLayersByName('test_flood')
layers = QgsProject.instance().mapLayers().values()
for point in iter_layer[0].getFeatures():
   counter += 1
   print(counter)
   for layer in layers:
       if layer.name() != iter_layer[0].name():
           for feature in layer.getFeatures():
               if (feature.geometry().contains(point.geometry())):
                       msgout = '%s, %s \n' % (point["pcds"], layer.name())
                       unicode_message = msgout.encode('utf-8')
                       output_file.write(unicode_message)
print('Finished')
output_file.close()

I am using QGIS 3.4.9-Madeira and my OS is MacOS 10.14.6 Mojave.


texts: XyTools QGIS plugin for attribute export

texts: I want to export an attribute table to Excel. I see XyTools is used for that (is there another plugin?). Since I installed QGIS from OSGeo4W, I have the required library. However, when I tell XyTools to export the table, I get the following Python error message: UnicodeEncodeError: 'ascii' codec can't encode character u'\xfc' in position 0: ordinal not in range(128) because my table contains nonASCII characters, which is not that unusual, especially for foreign languages. My question is: Is there a fix or workaround this issue or another plugin to export an attribute table to Excel?


texts: Getting statistics in cross sections at varying heights in LAS files?

texts:
I am looking to cross section my forest at varying heights. I have high return LAS data (8 points/m2 and average point spacing of 0.2m) and I am looking to do more than the simple 1, 2, 3 approach of vegetation classification within the LAS file.

I would like to have the flexibility to create multiple cross sections at specific height increments. Looking at, e.g., LASTools and ArcGIS tools, I can't do what I would like. So, I am looking at scripting something myself but I am running into personal knowledge roadblocks.

Firstly, when I create a LAS dataset in Arc each point that is created seems to represent 3500 LiDAR returns; when I run multi-point to single point this does not change.

Secondly, when I try and tile my problem using a fishnet approach the selected points for the current processing tile are not bounded by the tile and produce patterns of selected points well outside the tile (a multi-point property I'm sure).

How do I get at the data to perform my cross sections? Is there a way to play with LAS data outside of a LAS dataset?

UPDATE:

I have successfully sliced up my LAS data using some painful arcpy methods (see code). I am in the process now of running Jae's answer and the Fusion toolbox looks like the way to go in the future.

for file in os.listdir(in_las_folder):
    if file.endswith(".las"):
        print '...Working on LAS file {} of {}'.format(current, count)
        process = 2
        las = os.path.join(in_las_folder, file)
        las_name = os.path.splitext(file)[0]
        arcpy.CreateLasDataset_management(las, out_location_las + '\\ld_{}.lasd'.format(las_name), "NO_RECURSION", "", spatial_reference)
        lasd = out_location_las + '\\ld_{}.lasd'.format(las_name)
        arcpy.ChangeLasClassCodes_3d(lasd,"3 0;4 0;5 0")

        print '.....Creating Cross Sections'
        while process <= 38:
            print '........Working on heights between {}m and {}m'.format(process - 2, process)
            raster = out_location_height + '\\h{}_{}'.format(las_name,  process)
            if overwrite_height_model == 'NO':
                if arcpy.Exists(raster):
                    process += 2
                    continue
            print '.........Classify By Height'
            if process == 2:
                value = 20
                arcpy.ClassifyLasByHeight_3d(lasd, "GROUND", "19 1;20 2;21 4;22 6;23 8;24 10;25 12;26 14;27 16;28 18;29 20;30 22;31 24", "NONE", "false", "DEFAULT", "false")
            elif process == 26:
                arcpy.ChangeLasClassCodes_3d(lasd,"19 0;20 0;21 0;22 0;23 0;24 0;25 0;26 0;27 0;28 0;29 0;30 0;31 0")
                arcpy.ClassifyLasByHeight_3d(lasd, "GROUND", "19 24;20 26;21 28;22 30;23 32;24 34;25 36;26 50", "NONE", "false", "DEFAULT", "false")
                value = 20
            arcpy.MakeLasDatasetLayer_management(lasd, 'las_layer', [value])
            print '.......Creating raster'
            arcpy.LasDatasetToRaster_conversion('las_layer', raster, "ELEVATION", "BINNING MAXIMUM NONE", "FLOAT", "CELLSIZE", "0.5", "1")

            value += 1
            process += 2
        current += 1

        sys.exit()
print 'Creating full rasters'
process = 2
current = 1
while process <= 36:
    print '...Working on {} of {}'.format(current, count)
    for ras in arcpy.ListRasters():
        if int(ras.split('_')[1]) == process:
            height.append(ras)
    arcpy.MosaicToNewRaster_management(height, out_location_height, 'h{}'.format(process), spatial_reference, '32_BIT_FLOAT', '1', '1')
    for f in height:
        arcpy.DeleteFeature_management(f)
    process += 2
    current += 1


texts: Creating Vector Tile Package from basemap with ArcPy gives ERROR 001856: Cached scale doesn't match tiling scheme

texts:
I am trying to create a vector tile package from a basemap with arcpy, but get the using arcpy.management.CreateVectorTilePackage() causes me an error:

  ERROR 001856: Cached scale doesn't match tiling scheme.

My code is this:

import os

import arcpy

#set environment settings
arcpy.env.overwriteOutput = True

outputPath = "C://PROYECTOS//Pruebas_para_paquetes//Paquetes//"

minscale = 295828763

maxscale = 564

#   create a vector tile package,
#   using the same name as the map
aprx = arcpy.mp.ArcGISProject("C://PROYECTOS//Pruebas_para_paquetes//Datos_Prueba_Portal_ICA//Datos_Prueba//Mapas_ICA//Mapas_ICA.aprx")
map = aprx.listMaps()[3]

indices = r"C:\PROYECTOS\Pruebas_para_paquetes\Datos_Prueba_Portal_ICA\Datos_Prueba\Mapas_ICA\Mapas_ICA.gdb\Indices"

print("Packaging " + map.name)

arcpy.management.CreateVectorTilePackage(map,
                                        outputPath +'test.vtpk',
                                        "ONLINE",
                                        None,
                                        "INDEXED",
                                        minscale,
                                        maxscale,
                                        indices,
                                        "Mapa Base",
                                        "Mapa Base")


texts:
I have a similar question to:
How to join lines in QGIS for later using "lines to poly"?

I have a line shapefile that I wish to transform into a polygon shapefile.
After making some cleaning and closing some gaps, the lines have this look:

After that I've made "Vector>Geometry Tools>Lines to Polygons" but I've got this:

I've made a detection for dangles... and I found a few on the places that seems to have that kind of error. Also... the location of those points, seems to be on the places the "lines to poly" process fails.
I've tried to edit some of the dangles, snapping them to the nearest point to make sure there were no gaps... but running again the dangles detection, they're still there.
This is how all the shapefiles (plus the dangles) look like:

Is this easy to fix without "going to the grass"?
Why don't the dangles disappeared after snapping them to the correct place?

EDIT 1:
I was able to clean all the dangles... I was not happy with the snapping... and I've retried once more. I guess I was trying the wrong ones... I tried a very simple one... and one by one I was able to solve the more complex dangles, some were overshoot lines that were behind other segments and other were the intersections.
Well... after that, I was able to clean all the dangles!!!! :D

Now for the bad news... I got an error while making the Lines to Polygons:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:\PROGRA~1\QGISCH~1\apps\qgis\python\plugins\fTools\tools\doGeometry.py", line 348, in run
    success = self.lines_to_polygons()
  File "C:\PROGRA~1\QGISCH~1\apps\qgis\python\plugins\fTools\tools\doGeometry.py", line 519, in lines_to_polygons
    if inFeat.geometry().isMultipart():
AttributeError: 'NoneType' object has no attribute 'isMultipart'

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.4.0-Chugiak Chugiak, 8fdd08a
(...)]

Does this ring any bells to you guys?

(by the way... the first question was about the dangles... now it's about the lines to polygons process, however... the question is "Can't transform lines to polygons" should I close this and start another or may I continue here?)

EDIT 2:
user30184 asked for the Layer (no dangles at the moment):
Lines shapefile


texts: Nightly backup using Import XML Document arcpy ERROR 000791

texts:
I'm having some trouble using arcpy.ImportXMLWorkspaceDocument_management to import an XML document into a file geodatabase for a nightly backup routine. The code for the import is as follows:

 try:
     arcpy.ImportXMLWorkspaceDocument_management(backUP, backUpGDB, "DATA","DEFAULTS")

 except Exception as error:
     print error
 else:
     fgdbFin = datetime.now()
     print "Finished backup process at: " + str(fgdbFin)

The error I'm getting reads: "ERROR 000791: Cannot access local or remote geodatabase."

Oddly, if I try running the Import XML Workspace Document tool in Catalog, it runs without an issue. I'm running the script from a virtual desktop which only has 4GB of memory and it does seem to tie up most of my memory when running though I'm not sure why it would throw this error and not crash due to lack of resources.

Has anyone experienced this and if so are there any workarounds or suggestions for getting this to run in Python?


texts:
I am importing some very large rasters (in the tif format) into Postgres SQL (9.3) using raster2pgsql on the command line. Most of the rasters imported fine, some had to be converted to the img format and then imported fine, but some more have not imported at all. They appear to work, and the normal import line comes up, and then I get the following:

NOTICE:  Adding nodata value constraint
CONTEXT:  PL/pgSQL function addrasterconstraints(name,name,name,boolean,boolean,
boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean)
line 53 at RETURN
NOTICE:  Adding out-of-database constraint
CONTEXT:  PL/pgSQL function addrasterconstraints(name,name,name,boolean,boolean,
boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean)
line 53 at RETURN
NOTICE:  Adding maximum extent constraint
CONTEXT:  PL/pgSQL function addrasterconstraints(name,name,name,boolean,boolean,
boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean)
line 53 at RETURN

and the import appears to make no further progress.

I have searched the metadata and found no differences at all from the files that worked, the same command was used to import them all:


texts: Polygon to Raster Cells Are Stretched Out

texts: I am using ArcGIS Pro, polygon to raster tool. I am converting a large set of polygons (over 100,000) into a raster, using a cell size of 0.0016 (the default is 0.0098). When the operation is complete, it creates cells that are very stretched out horizontally. The cells seem to be the correct vertical distance, but are way too wide. I tried the operation with the same inputs but on a small selection of polygons, and it worked fine.


texts: Can't Get my legend on a swipe layer map

texts:
I'm having trouble getting a legend on my swipe layer map. I'm new to Javascript and took most of the code from the API reference but can't get the legend to show up for either layers. I'm Hoping someone can take a look at what i'm doing and give me some advice what i'm doing wrong.

  <!DOCTYPE HTML>
<html>
<head>
  <meta charset="utf-8">
  <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
  <title>Layer Swipe</title>
   <link rel="stylesheet" type="text/css" href="http://js.arcgis.com/3.11/esri/css/esri.css">
<style>
html, body, #map,.map.container {
  padding:0;
  margin:0;
  height:100%;
  width:95%;
}
#legendDiv{
background-color: #fff;
position: absolute !important;
z-index: 99;
top:01%;
right:0%;
 }
 </style>
 <script src="//js.arcgis.com/3.11/"></script>
 <script>

require([
  "esri/map",
  "esri/dijit/LayerSwipe",
  "esri/arcgis/utils",
  "dojo/_base/array",
  "esri/dijit/Legend",
  "dojo/domReady!",
], function(
  Map, LayerSwipe, arcgisUtils, array, Legend
)  {

  var mapDeferred = arcgisUtils.createMap("e80b5668e7ea478c986dff4bfbe7112a", "map");

  mapDeferred.then(function(response){

    var id;
    var map = response.map;
    var title = "Inundation Post Dike Placement";

    //loop through all the operational layers in the web map
    //to find the one with the specified title;
    var layers = response.itemInfo.itemData.operationalLayers;
    array.some(layers, function(layer){
                  id = layer.id;

    });
    //get the layer from the map using the id and set it as the swipe layer.
    if(id){
      var swipeLayer = map.getLayer(id);
      var swipeWidget = new LayerSwipe({
        type: "vertical",  //Try switching to "scope" or "horizontal"
        map: map,
        layers: [swipeLayer]
      }, "swipeDiv");
      swipeWidget.startup();
    }
  });

  var legend = new Legend({
      map: map,
      layerInfos:(arcgisUtils.getLegendLayers(response))
  }, "legendDiv");
  legend.startup();
});


texts: Take parallel lines and automatically move them apart

texts:
Rather than collapsing parallel lines, I am looking for a process that would separate them.

There are a number of lines on the network that run too close together to distinguish as separate features and I would like to 'push' them apart from each other.

The image below has an inset that shows the features where they visually merge, but the full picture shows the lines at a larger scale.

Anything?


texts:
i use two seperate datasource for my layers:

sql database for main map data, this is using EPSG:900913
osm file for some extra features, this data use EPSG:4326 projection

tiles for both layers are created on server side.
what i need to do to generate tiles for EPSG:900913 projection from data in EPSG:4326 projection?
is it somewhere in tilecache configuration? in definition map or layer in osm file, or somewhere else?

tilecache.cfg

[test]
type=Mapnik
mapfile=/var/www/tc/osm/test.xml
tms_type=google
layers=global

test.xml

<Map background-color="#00000000" srs="+proj=latlong +ellps=WGS84 +datum=WGS84 +no_defs" >
  <Style name="style">
    <Rule>
      <PolygonSymbolizer fill="#ff0000" />
      <LineSymbolizer stroke="#00ff00" stroke-width="100" />
    </Rule>
  </Style>
  <Layer name="global" srs="+proj=latlong +ellps=WGS84 +datum=WGS84 +no_defs" status="on">
  <StyleName>style</StyleName>
  <Datasource>
      <Parameter name="type">osm</Parameter>
      <Parameter name="file">/var/www/tc/osm/test.osm</Parameter>
    </Datasource>
  </Layer>
</Map>

test.osm

<osm version='0.6' upload='true' generator='JOSM'>
  <node id='-186' action='modify' visible='true' lat='50.0' lon='50.0' />
  <node id='-184' action='modify' visible='true' lat='50.0' lon='-50.0' />
  <node id='-182' action='modify' visible='true' lat='-50.0' lon='-50.0' />
  <node id='-180' action='modify' visible='true' lat='-50.0' lon='50.0' />
  <way id='-91' action='modify' visible='true'>
    <nd ref='-186' />
    <nd ref='-184' />
    <nd ref='-182' />
    <nd ref='-180' />
    <nd ref='-186' />
  </way>
</osm>

map.html

...
var   map = new OpenLayers.Map ("map", {
                controls:[
                    new OpenLayers.Control.MousePosition(),
              numZoomLevels:zoomLev,
                maxResolution: 156543.0339,
                units: 'm',
                projection: new OpenLayers.Projection("EPSG:900913"),
                maxExtent: new OpenLayers.Bounds(-20037508.34, -20037508.34, 20037508.34, 20037508.34),
                displayProjection: new OpenLayers.Projection("EPSG:4326")
            } );

var test = new OpenLayers.Layer.WMS( "test", "http://localhost/tc/tilecache.cgi?",
    {layers: 'level_1', format: 'image/png'}, { isBaseLayer: false} );
...


texts: How to use shortcuts in dockable windows?

texts:
I have developed a custom dockable window with textboxes, comboboxes and buttons for use in ArcMap. After the user clicks on a button in the dockable window a function is invoked.
How can I invoke the function/button_click-event using short cuts. Where do I have to wire the event?:

this.KeyDown +=new KeyEventHandler(MyDockableWindow_KeyDown);

Where do I have to impement the event?

private void MyDockableWindow_KeyDown(object sender, KeyEventArgs e)
    {
        switch (e.KeyCode)
        {
            case Keys.F1:
                button1.PerformClick();
                break;
            case Keys.F2:
                button2.PerformClick();
                break;
            case Keys.F3:
                button3.PerformClick();
                break;
            case Keys.F4:
                button4.PerformClick();
                break;
        }
    }


texts:
I've checked the other answers here and none of the answers worked for me...

It's a "simple" matter of finding a layer by name in a map and changing it's datasource.

Here's my script:

    #-- New Source
    newSource = os.path.join(currDir, lyrZones + "_Dep" + dirNum + ".shp")

    #-- Validation
    if arcpy.Exists(newSource):
        print "The new source exists!"
    else:
        print "Not found!"

    #-- Find the layer and change its source
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.name == lyrZones:
            print "Layer found!"
            print layers.dataSource
            print "  Replacing source..."
            lyr.replaceDataSource(currDir, "SHAPEFILE_WORKSPACE", newSource)
            print lyr.dataSource

I run it from a python script outside ArcGIS.

You can see that the source exists, the layer is found, I can access it's datasource by printing it, but the lyr.replaceDataSource causes an error (Unexpected error... doesn't help much...)

I've been working at it for 4-5 hours now and I'm short of solutions to try. Any idea what can cause the crash?


texts: Merging polygons based on separate shapefile touching each polygon using ArcGIS Desktop?

texts:
In the image provided, I have drawn few polygons and a shapefile in red colour above them touching those polygons. Is there any way by which we could merge those polygons which, the red colour shapefile is touching.

At the end what I am expecting is that suppose the red polygon is touching 3 polygons, they should make it a single polygon(just the polygons which are touched by red polygon/shapefile)


texts:
I have an image that I want to thin, after calling the image on the grass and executing the r.thin command, it was not thinned and the output image is as follows:
image original:

And output image of the command r.thin:

Why do you think thinning is not done and what should I do?


texts: Iterating Intersect_analysis over polygon shapefile using ArcPy?

texts:
I am using ArcGIS 10.3.1 for Desktop and I have two shapefiles:

A set of stream (vector line) data
A fishnet of polygons

I want to know, for each polygon in that fishnet, what the total length of streams is within that polygon. So far, I've tried iterating the intersect_analysis tool over the fishnet table, summing the lengths of the output and writing that to the relevant field, but each iteration performs the intersect on the entire grid, not the selected polygon in that grid.

Code:

import arcpy

#set some environmental variables
arcpy.env.workspace = "C:\Users\Alex\Documents\SoilMoisture"
arcpy.env.overwriteOutput = 1
#define variables for all necessary feature classes and fields

streams = "pathto/str_split_lamb.shp"  ## shapefile of stream network split at vertices, in Lambert projection, with Length value field pre-calculated
rainGrid = "pathto/net_test.shp" ## polygons for rain cells, with empty "DRN_LENGTH" field

cursor = arcpy.UpdateCursor(rainGrid) ## will iterate across every row of the rainfall grid feature

for row in cursor:

    intStreams = arcpy.Intersect_analysis([streams,rainGrid],"output.shp","ALL","","LINE") ## intStreams becomes all stream segments inside polygon

    strCursor = arcpy.SearchCursor(intStreams) ##iterate over streams in intersection

    totalLen = 0
    for rowmore in strCursor:
        totalLen += rowmore.length ## sum all length values into one variable

    print totalLen # this is just here to check what's going on

    row.DRN_LENGTH = totalLen ## add value to field

    cursor.updateRow(row)

print "Done!"

Any ideas?


texts: Locate roads that are close to steep slopes which risks landslides

texts:
I've a query, I'm trying to locate all the roads which are basically at risk of erosion (landslide?) when it's raining heavily for a municipality.

What I have to work with is this:

A shape layer containing lines for the roads (centerlines)
Contourlines for the terrain
Raster with elevation model
ArcGIS Advanced Desktop with 3d and spatial analyst

My first thought was to create some lines perpendicular to the roadlines, in a range of for instance 20 meters, and there after count how many contourlines i cross. This is pretty lowtech, but so far it's the best I've come up with.

Does anybody have a smarter approach? Or perhaps tried doing something similar?


texts:
I want to make a Data Driven export to PDF of the layout using arcpy, and later add to that PDF extra pages from PDFs from my local computer.

My problem is that I don't actually want to make the exports of the PDFs on my computer using the arcpy.mapping.ExportToPDF , I would only want to store that page (the data driven layout) in the memory, and later on, add the extra pages.

I use this code:

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, '')[0]

for pageNum in range(1, mxd.dataDrivenPages.pageCount+1):
    pageName = mxd.dataDrivenPages.pageRow.getValue(nameField)
    arcpy.mapping.ExportToPDF(mxd, r"C:\TEST\tempExport{}.pdf".format(pageNum)) # This is the temporary export that I want to be stored in memory, not on my local machine
    pdfDoc = arcpy.mapping.PDFDocumentCreate(r"C:\TEST\FinalPDF{}.pdf".format(pageNum))
    pdfDoc.appendPages(r"C:\TEST\tempExport{}.pdf".format(pageNum))
    pdfDoc.appendPages(r"C:\TEST\local{}.pdf".format(pageNum))
    pdfDoc.saveAndClose()


texts: r.report giving only per-class cell count

texts:
I'm trying to extract the count of each pixel value in a raster, like in this question: How to extract pixel value counts from a raster in QGIS?
However, when I run r.report (from QGis (2.10) Processing toolbox) I get a result like this:

Even if my settings in r.report are the same as the ones in that question i linked earlier:

I tried to change the settings but I can't get the result I want. I don't know if I'm missing something obvious, but I've been looking around for a while and found no explanation for this.


texts: Duplicating multiple shapefiles in multiple locations

texts:
I have created a design for a well pad with two well heads and an access road that I would like to duplicate in multiple (20) locations (i.e., the same site design will be used at multiple sites.) Because the design involves points (well heads), lines (roads) and polygons (well pads), where the spatial arrangement is important, I have multiple feature classes and I don't want to simply copy one, drag to a new location, then copy another and try to fit it back where it should go. Instead, I would like to copy the entire design, consisting of three feature classes--one of points, one of lines, and one of polygons--and paste it in a new location. Provided I am able to drag all of the elements together, I can use aerial imagery to line up the features where they should go on each site.

I'm indifferent as to whether I have multiple copies of each feature in the same feature class (i.e., 20 copies of the well pad) or multiple feature classes (i.e., 20 feature classes with one well pad).

Any suggestions?

Using ArcGIS 10.3, Desktop, Basic license, with Spatial Analyst Extension.


texts: Batch export of layers to PDF producing blank PDFs

texts:
Based on Exporting set of layers into multiple PDF (or Adobe Illustrator) files, one at a time? I was able to alter and execute code that is trying to export 30+ layers from an MXD I have.

import arcpy

def turn_off_layers(mxd, df):
    for lyr in arcpy.mapping.ListLayers(mxd, "", df):
         lyr.visible = False
    mxd.save()

def print_layers(mxd,df):
    for i,lyr in enumerate(arcpy.mapping.ListLayers(mxd, "", df)):
         lyr.visible = True
         mxd.save()
         arcpy.mapping.ExportToPDF(mxd, r"Y:\Desktop\Google Drive\DTI\Figures\%s.pdf" % lyr.name)
         lyr.visible = False
         mxd.save()

if __name__ == "__main__":
    mxd = arcpy.mapping.MapDocument(r"Y:\Desktop\Google Drive\DTI\DTI_Onondaga_100915.mxd")
    df = arcpy.mapping.ListDataFrames(mxd, "")[0]
    turn_off_layers(mxd, df)
    print_layers(mxd,df)

Unfortunately, the code runs without error and names the PDFs correctly, but every PDF is blank.  Can anyone see where I am going wrong?
I am using ArcMap 10.2.2 on a 2011 iMac, running Parallels 9 in Cohesion view.

EDIT:  Edited code based on Farid's and PolyGeo's suggestions:

import arcpy

def turn_off_layers(mxd, df):
    for lyr in arcpy.mapping.ListLayers(mxd, "", df):
         lyr.visible = False

def print_layers(mxd,df):
    for i,lyr in enumerate(arcpy.mapping.ListLayers(mxd, "", df)):
         print lyr.name
         lyr.visible = True
         df.extent = lyr.getExtent()
         arcpy.RefreshTOC()
         arcpy.RefreshActiveView()
         arcpy.mapping.ExportToPDF(mxd, r"Y:\Desktop\Google Drive\DTI\Figures\%s.pdf" % lyr.name)
         lyr.visible = False

if __name__ == "__main__":
    mxd = arcpy.mapping.MapDocument(r"CURRENT")
    df = arcpy.mapping.ListDataFrames(mxd, "")[0]
    turn_off_layers(mxd, df)
    print_layers(mxd,df)

AND THE EDITED CODE WORKS!


texts: Extracting raster values to point features with reduce regions. Error: “User memory limit exceeded”

texts:
I am very new to Earth Engine and Javascript.

I'm trying to get the cumulative cost distance between some areas of interest and some tide gauges in the coastal U.S. To do this, I first calculated the pixel cost based on a national elevation map. Then, I calculated the cumulative cost using to the built-in function cumulativeCost. That went all pretty well. Now I'm trying to extract the cumulative cost values at the position of the tide gauges. To do this, someone suggested me to use the reduceRegions method. I've tried the following code but unsuccessfully.

I'm going to post my whole code so that is replicable. Please note that the part with which I have a problem is the second one.

// IMPORTS
var imageCollection =
ee.ImageCollection("users/brazzolanicoletta/sourceRasters"),
sourceVis = {"opacity":1,"bands":["b1"],"min":1,"max":1,"gamma":1},
DEMVis = {"opacity":1,"bands":
["elevation"],"min":-73.50744474674659,"max":374.555654458347,"gamma":1},
imageVisParam = {"opacity":1,"bands":
["elevation"],"min":-0.02534169006347656,"max":3.6601884765625,"palette":
["0345ff","000000"]},
imageVisParam2 = {"opacity":1,"bands":
["elevation"],"min":-0.02534169006347656,"max":3.6601884765625,"palette":
["0345ff","000000"]},
cosVisParam = {"opacity":1,"bands":
["cumulative_cost"],"max":4170.014060708561,"palette":
["ff0303","efff05","4eff05","002bff","ff01f7","000000"]},
imageVisParam3 = {"opacity":1,"bands":
["cumulative_cost"],"max":4028.1446098656247,"gamma":1};

//get IDs for images in image collection
var getID = function(image){ return image.set('ID', image.id());};
var okID = imageCollection.map(function(image) { return image.set('ID',
image.id());});

// Set general estethic parameters
var dem_vis = {bands:"elevation", min:0, max:0.05,
palette:"#0345ff,#000000"};
var cost_vis = {bands:"cumulative_cost", min:0, max:10000,
palette:"ff0303,efff05,4eff05,002bff,ff01f7,000000"}

// PART 1: Cumulative Cost based on source rasters
 //import elevation map
var dem = ee.Image('USGS/NED');

// pixel cost calculation
var elThreshold = ee.Number(5); //set elevation threshold
var subDEM = dem.updateMask(dem.lt(elThreshold)); //mask pixel above
elevation threshold
var costDEM = (subDEM.add(30)).divide(1000); //calculate the cost of each
pixel (height + width pixel (30m)) in km

// Add DEM to the map
Map.addLayer(costDEM, dem_vis, "SRTM");

// Cumulative cost
var calcCumCost = function(img) {
return costDEM.cumulativeCost({
source:img,
maxDistance:1E5});
}; //write a function that perform the cumulative cost calculation for each
image given the cost of the pixel

var demCost = ee.ImageCollection(okID.map(calcCumCost)); // caulcuate
cumulative cost for each source raster in the image collection

// PART 2 - Reduce Region: extract cumulative cost for tide gauges

var tideGauges =
ee.FeatureCollection('ft:1e1ik7ZklKbRSRVS50Ml_prHBTZ0WbNgW73fw7Ald');
//import fusion table of tide gauges

// WORK IN PROGRESS
// Empty Collection to fill
var ft = ee.FeatureCollection(ee.List([]));

// function that extract values from cumulative cost rasters and reduce it for points region
var fill = function(img) {
 // gets the values for the points in the current img
var ft2 = img.reduceRegions(tideGauges, ee.Reducer.first(),30);
// set ID
var ID = ee.Feature(null, {'id':img.id()});

// writes the ID in each feature
 var ft3 = ft2.map(function(f){return f.set("id", ID)});
// merges the FeatureCollections
return ft.merge(ft3);
};

// Apply the function to each image in the ImageCollection
var newft = ee.FeatureCollection(demCost.map(fill));
print(newft, 'Potentially: region-reduced cost');


texts: SLD for Line-Features with gap at start and end

texts:
I'm trying to define an SLD for Line-Features (Using for WMS in Geoserver). The signature of the line mustn't go from startpoint to endpoint, but has to start 5 pixels AFTER the startpoint and has to end 5 pixels before the endpoint.

O______O

Background:
I want to visualize the start-/endpoint from lines with a ring-symbol. The problem is, that the line mustn't intersect the ring-symbol. Normally, I would fill the ring-symbol, nevertheless, underlying polygon-features shouldn't be covered by the filling of the ring.