texts: Finding upstream and downstream villages along river?

texts:
I have a DEM extract of river and village centroid. How do i create an attribute table that shows which village is upstream of other or downstream to another.

Can someone guide me the correct tools in arcgis?

Also im not quit comfortable with programming.


texts:
Given a file my_file.shp such those from Gadm e.g. US administrative divisions, where the coordinate reference system is latitude/longitude with WGS84 datum.

Input could actually be:

Shapefile,
ESRI personal geodatabase,
ESRI file geodatabase,
Google Earth .kmz,
R (SpatialPolygonsDataFrame)

How to convert this input into a WGS84 PostGIS dump ?


texts: How to intersect 2 polygon layers with a given tolerance for features in different layers only?

texts:
I can intersect 2 layers with polygons using ArcMap. However I can specify only overall tolerance no matter what layer a feature is in (even same layer). So I either end up with lots of slivers if I don't tolerate snapping, or I end up collapsing small features from the first layer that I'd want to keep.

The problem is that I would like for features let's say less than 10 m apart in the first layer to stay as they are if they happen to be contained within some polygon of the second layer. However I'd like for edges of polygons of the 2nd layer to snap to those in the first layer if they are within 10 m.

I'm open to a solution using any tool be it ArcMap or some python library/script.

Here is the rough idea of what I want.


texts: Copying Features from One ArcMap window to another, without losing Related Tables

texts:
In our current workflow we are working in a Personal Geodatabase with few features and related tables , many persons will be working in different schemes.

When we need to merge these schemes\Pgdb to a single output , we are copying feature class to feature class . Obviously this changes the Object ID in the resulting output Geodatabase and we lose the related table Information.

Is there any way to copy feature classes from different inputs Pgdbs and still preserve the related table ?


texts: How to display the join field in a QGIS join?

texts:
Upon joining tables in QGIS (2.8.2) I am unable to view the join field in the resulting join.  Additionally, the join field is not listed in the target table's Layer Properties, even though the other join table fields are listed there.

I recognize that the target field is displayed, but there are times when I would like to also see the join field.  How do I do this?


texts:
I'd like to do zonal statistics / terrain analysis with polygon overlays - I achieved to get min / max like so:

library(raster)
library(rgdal)

alt <- getData('alt', country = "AT")
gadm <- getData('GADM', country = "AT", level = 2)
gadm_sub <- gadm[1:3, ]

plot(alt)
plot(gadm_sub, add=T)

extract(alt, gadm_sub, fun = max, na.rm = T, small = T, df = T)

But for slope / aspect I get this error which I don't grasp how to prevent:

extract(alt, gadm_sub, fun = function(x) terrain(x, opt = "aspect"), na.rm = T, small = T, df = T)
# > Fehler in FUN(X[[1L]], ...) : unbenutzte(s) Argument(e) (na.rm = TRUE)


texts: Seeking training resources in Python as applied to geodatabase management?

texts: I've taken ESRI database courses, which have been helpful to my work but I also need training in Python as applied to geodatabase management?


texts:
In a Arcmap raster weighted overlay, if an input data set is designated as "restricted" will overlay calculations not be performed on the "restricted" areas, and those areas will be included in the output? Or will The overlay function ignore the "restricted" data altogether?

I guess I'm just trying to determine how "restricted" could be compared with "NoData". I have heard "restricted" is included in the output, and I have heard it is ignored altogether, and thus want to understand which is more accurate.


texts: Using IF/ELSE Statement for field calculation in ArcPy?

texts:
I would like to calculate the slope value in %, but I don’t want to have any negative values. If there are negative values in the slope field, I would like to multiple them by -1 and if not, multiple them by +1 (or do nothing).

My input is a line shapefile. First I added a field called "slope", then I calculate the slope based on other fields.

arcpy.AddField_management(ShapeFileLine, "Slope", "DOUBLE", 10, 3)
arcpy.CalculateField_management(ShapeFileLine, "Slope", '((((!END_Z! - !START_Z!)/!LENGTH_3D!))*100)!',"PYTHON")

What I should do if I want I only have positive values? In this case they are all negatives, but sometimes the positive and negative values are mixed


texts: iPhone app to record 3D location information

texts:
Does anyone know of an app that would run on an iPhone or iTouch that would record the X,Y and Z information from an external GPS?

I am using a Dual XGPS150 Bluetooth receiver.  The GPS will record X,Y,Z info, but I am having issues with an app that will record the info, and provide it back to me in a useful format.  I have tried GISRoam, but it will not record the Z information.  I have looked at other apps, but I am not sure if it will record the Z info and allow me to export it into a format that can be used in a GIS solution like QGIS.  Does anyone have any experience with the WolfGIS, Map My Run, or any of the other IOS apps?


texts: How to automaticaly draw polygons inside polygon?

texts:
I would like to know if it is possible to automatically draw polygons patterns in a polygon.
With more details I would like to go from this :

To this :

The large rectangles are 6.90m*31.05m and can be reduce to 6.90m*1.15m.
The spacing between each array is 11.5m, ideally I would like to be able to adapt this value.
I understand how to create the border, but not how to generate these polygons inside the first one.
I think I am down to use the "model" and the PyQGIS parts of QGIS but I don't know where to begin.
I think I also have an optimisation issue about this problem.


texts:
I have a large set of .csv-files containing geographic information. Every dataset consists out of 2 .csv-files! This is my folder-structure. All files in Tracks are formatted equally and all in Routes too, obviously.

------- Tracks
-- 1.csv
-- 2.csv
-- 3.csv
-- ...
------- Routes
-- 1.csv
-- 2.csv
-- 3.csv
-- ...

I can import them all for themselves, but that would take days of time (2000 datasets ...).

So what I would like to do is to load all those .csv-files at once and group them. So after the import I would like to have a lot of layer groups named 1,2,3 (like the files before the comma) and every group contains the layers representing the two files. Like this:

---- Group 1
-- tracks1
-- routes1
---- Group 2
-- tracks2
-- routes 2
---- Group 3
-- ...

Is that possible in any way? I cannot write python code, so I'm hoping for a plugin ;)

Thanks a lot!


texts: SLD filter with embedded strTrim function in Geoserver

texts:
Using Geoserver 2.2.2 (GeoTools Version 8.4) on Windows 7, I'm trying to create a SLD that trims the leading and trailing spaces of a string property (using the strTrim function) and compares it to a given string pattern (using PropertyIsLike) for filtering. Here is the code:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
xmlns="http://www.opengis.net/sld"
xmlns:ogc="http://www.opengis.net/ogc"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">

  <NamedLayer>
    <Name>style1</Name>
    <UserStyle>
      <Title>style1</Title>
      <Abstract>style1</Abstract>
      <FeatureTypeStyle>
        <Rule>
          <Name>Fauche</Name>
          <Title>Fauche</Title>
          <ogc:Filter>
            <ogc:PropertyIsLike wildCard="*" singleChar="#" escape="!">
              <ogc:Function name="strTrim">
              <ogc:PropertyName>type</ogc:PropertyName>
              </ogc:Function>
              <ogc:Literal>Testword</ogc:Literal>
            </ogc:PropertyIsLike>
          </ogc:Filter>
          <PolygonSymbolizer>
            <Fill>
              <CssParameter name="fill">#FF8C00</CssParameter>
              <CssParameter name="fill-opacity">0.6</CssParameter>
            </Fill>
          </PolygonSymbolizer>
        </Rule>

  </FeatureTypeStyle>
</UserStyle>

When validating the SLD, I get the following error:

  cvc-complex-type.2.4.a: Invalid content was found starting with
  element 'ogc:Function'. One of
  '{"http://www.opengis.net/ogc":PropertyName}' is expected.

Where is the mistake?


texts:
The resulting shapefile from my ArcGIS model is just disappearing after the model is finished running.

I see it for a short moment in catalogue but it doesn't show up in table of contents (I turned on add to display in the model).

Then if I refresh the folder that contains it in the catalogue, it also disappears from there.


texts:
I am trying to make a  systematic cluster sampling pattern for data collection. The sampling consist of the points placed over the landscape in 6 km span (red). From this major point, I need to create 3x3 sampling spots in 200 m span (grey), in direction towards the east and south, creating L shaped cluster of points.

I was wondering, how to make this sampling scheme in ArcGIS using ArcPy?

It is pretty straightforward to generate a set of equally distributed points over the landscape using Create Fishnet tool and set Create Label Points as TRUE. Cell height and cell width will simply be 6 km.

But do you have any suggestion, how to create set of grey points located towards east and south, in 200 m distance from each other?

I honestly even don't know where to start...

My sample code to generate red grid points:

import arcpy, os
from arcpy import env

# Allow files to overwrite
arcpy.env.overwriteOutput = True

# Set coordinate system of the output fishnet
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(3067)

# Define working directory
inWD = "C:/Users/data/raw"

# Define variables:
inFC = os.path.join(inWD, "myFC.shp")
outFC = os.path.join(inWD, "myFC_grid.shp")

# Create fishnet
cellSizeWidth = 6000 # 6 km
cellSizeHeight = cellSizeWidth
labels = "LABELS" # Create points in the middle of squares
desc = arcpy.Describe(inFC)  # get describe object to define extent coordinates

# Process: Create fishnet
arcpy.CreateFishnet_management(out_feature_class = outFC,
                               origin_coord = str(desc.extent.lowerLeft),
                               y_axis_coord = str(desc.extent.XMin)  + " " + str(desc.extent.YMax + 100),
                               cell_width = cellSizeWidth,
                               cell_height = cellSizeHeight,
                               number_rows = "0",
                               number_columns ="0",
                               corner_coord = str(desc.extent.upperRight),
                               labels = "LABELS",
                               template = "#",
                               geometry_type = "POLYLINE")

print "basic fishnet created"

# Import newly generated point grid data ""
red_points = os.path.join(inWD, "myFC_grid_label.shp")

# Make new freature class to add new points
grey_points = os.path.join(inWD, "Lcluster.shp")

How to add clusters points to grey dataset?


texts: Problems with features on ImageStatic source

texts:
I am trying to draw a grid on top of an ImageStatic source in OpenLayers but I am getting extra grid squares on the right side of my image.  Why do rows 2 & 3 have 5 columns instead of 4?
https://jsfiddle.net/mortac8/qh83vdp0/68/
var extent = [0, 0, 390, 260];
var projection = new ol.proj.Projection({
  code: 'test-image',
  units: 'pixels',
  extent: extent
});

var imageStaticSource = new ol.source.ImageStatic({
  url: "",
  projection: projection,
  imageLoadFunction: imageLoadFunction,
  imageExtent: extent
});

var source = new ol.source.Vector();

let topPx = 260;
let rightPx = 390;
let left = 0;
let size = 128;
let right = left + size;
let counter = 1;
while(topPx > 0) {
  while(left < rightPx) {
    console.log("adding box " + counter++);
    var ring = [[left, topPx], [right, topPx], [right, topPx-size], [left, topPx-size], [left, topPx]];
    var geometry = new ol.geom.Polygon([ring]);
    var feature = new ol.Feature({
      geometry: geometry
    });
    source.addFeature(feature);
    left = left + size;
    right = left + size;
  }
  console.log("NEW ROW");
  topPx = topPx - size;
  left = 0;
}

var vectorLayer = new ol.layer.Vector({
  source: source
});

var map = new ol.Map({
  layers: [
    new ol.layer.Image({
      source: imageStaticSource
    }),
    vectorLayer
  ],
  target: 'map',
  view: new ol.View({
    projection: projection,
    center: ol.extent.getCenter(extent),
    zoom: 0,
    maxZoom: 3,
    minZoom: 0
  })
});


texts: Open Layers 4 fit to extent zoom too much out

texts:
Working with OpenLayers 4.2.0, I have the extent [-45930163.095729224, -3380418.677804029, -5390354.729126186, -3011487.2342575104] which is a box around Santa Catarina state, Brazil.

But instead of having the map center in this box, I'm actually being centered in Australia, which would be an outer extend, if that is the right word for it, such as the map center as:

My code is (ignore the extent visibility and the condition for "Infinity"):

        extent = [-45930163.095729224, -3380418.677804029, -5390354.729126186, -3011487.2342575104];
        if (extent[0] != "Infinity") {
            map.getView().fit(extent, {size: map.getSize(), padding: [15, 0, 0, 0]});
        }

I've already tried to use a size different than map.getSize(), but all I get is the map zoomed in Australia.

Since I can't comment, @bartvde: It didn't helped too much. Kind of zoomed in, but the center is still in Australia.


texts: Why I'm getting "Runtime error : Object: CreateObject Layer invalid data source" while adding shapefile using ArcPy?

texts:
I am trying to add shapefile using ArcPy, but each time my script is ending in showing runtime error? I am using ArcGIS 10.

Here is my script

import arcpy
from arcpy import env

Master="G:\path\To\ShapeFile.shp"

mxd = arcpy.mapping.MapDocument("CURRENT")

df = arcpy.mapping.ListDataFrames(mxd,"*")[0]
layer = arcpy.mapping.Layer(Master)
arcpy.mapping.AddLayer(df, layer,"BOTTOM")

arcpy.RefreshActiveView()
arcpy.RefreshTOC()

This script works fine when I add .lyr files, but why it isn't working for .shp


texts: UnboundLocalError: local variable referenced before assignment

texts:
I have following simple function to get percent values for different cover types from a raster. It gives me following error:
UnboundLocalError: local variable 'a' referenced before assignment

which isn't clear to me. Any suggestions?

def percent(ras):
    temp = arcpy.SearchCursor(ras,"","","Value;Count","")
    for row in temp:
        theval = row.getValue("Value")
        thecount = row.getValue("Count")
        if row.getValue("Value") == 1:
            a = thecount
        if row.getValue("Value") == 2:
            b = thecount
        if row.getValue("Value") == 3:
            c = thecount
        if row.getValue("Value") == 4:
            d = thecount
    return a,b,c,d


texts: Google Earth Engine: how to negate the .and() operator on an image with NULL values

texts:
How does one negate a ".and()" operation on an image in GEE if this image does not contain 0/1 values but NULL/1 values?

From the documentation:

  not() - Returns 0 if the input is non-zero, and 1 otherwise. Arguments:
  this:value (Image):
      The image to which the operation is applied. Returns: Image

For example if I want to retreive all loss in the Hansen tree loss dataset that does not overlap with a tree cover density >= 30. I use this code:

var gfw = ee.Image('UMD/hansen/global_forest_change_2018_v1_6');
var tc = gfw.select('treecover2000');
var loss = gfw.select('loss');
var lossyear = gfw.select('lossyear');
var frst = tc.gte(30);
var not_overlapping_loss = lossyear.and(frst.not())

which works fine. But when I try the same code replacing "frst" by an "asset" I uploaded that contain only "NULL/1" values, it seems the resulting image contains only NULL values as the (not existing) 0's are never converted to 1's...


texts:
I converted the file I exported from OSM to PostGIS using osm2pgsql. I have opened the table's outputs to QGIS and I saw data are not complete compare to OSM editor view. Most of the roads are broken and out of align.

How do I complete them?


texts:
I am using QGIS 2.4.0 Chugiak version.

Is there an option to save selected features from a shapefile?

I couldn't find such an option in this version. I am using stand alone application.


texts: Split by Attributes tool exporting to DBF files only?

texts:
I am using the "Split by Attributes" tool in ArcGIS 10.6 to go through a shapefile, and split it up into new files by an attribute value. The tool runs perfectly, however I notice it only exports to DBF files. Shapefiles or CSVs would be preferable, but I don't see any option to specify the file type.

Can this be changed and, if not, is there a way to split up a shapefile by attribute value, and export it to .shp or .csv files?


texts: Union with No_Gaps puts "0" into blank fields using ArcPy?

texts:
wanted to overlay two file geodatabase-based featureclasses where one has a smaller extent than the other, but preserve all input features in the output. I used arcpy.Union_analysis. In order to maintain feature geometry present in even only one of the two input FCs, I use the "NO_GAPS" option. The documentation seems (to me) to state that such features will have the full/combined attributes of both feature classes, but that values will be "blank" (actually says: This feature will have blank attributes.)

Problem is that the values for numerical attributes of these features are not actually "NULL" in the output. They are zero.

Is this correct or is there any way to ensure that "NULL" values are produced?

Zero is a legitimate value, so having that in the output is a failure for my use case.


texts: What is the best way to perform a diff operation on Postgis features/tables or Geotools objects? Are there any automatic methods I can use? If not, what would be the optimal way to implement your own diff? For example I have two versions of a Shapefile and I import them into two separate PostiGIS tables, what algorithm should I consider now - how to compare a geometry column and get a diff result? I would be grateful for any advices or links to additional sources/articles concerning this matter.


texts: Custom contour color gradient using RGB values from QGIS

texts:
I'm using QGIS 3.4.1, and downloaded the SRTM data for the area of my interest, and added the TIFF format SRTM data file to my project using "add raster layer" option.

Here, I wish to change the contour color gradient for SRTM data from property-symbology menu. More precisely, I wish to color the height in terms of RGB value. For the project of mine, I only need to use red and green colors.

Following is the detail:

Formula is (Red Channel) + (Green Channel * 255).

For example: if a pixel value is RGB (171,1,0) then the height is (171) + (1 * 255) = 426

Another example: if a pixel value is RGB (107,2,0) then the height is (107) + (2 * 255) = 617

How can I differ my contour color for my raster layer (SRTM TIFF file) in the same way with those example?


texts:
I'm looking for a free download of a shapefile about forest roads in Europe (I need forest roads of Austria and Slovenia).

Here in Italy you can free download these kind of shp, but actually I can't find a link or a website for downloading austrian and slovenian forest roads, probably because many websites are written in their own languages which I don't understand (fortunately I found highways and urban roads shapefiles, and now I miss only forest roads).

Any suggestion?


texts: Exporting feature colors from categorized polygon layer to Excel diagram using QGIS?

texts: I have an xxx.shp layer with polygones which is categorized and colored according to the names of management companies. I need to make a circle diagram for this layer in Excel. Is it possible to make this by any plugin?


texts: How can I specify the CRS of a raster layer in pyqgis?

texts:
I manipulating netcdf raster data from the QGIS python console.
The metadata of my netcdf file does not provide the Reference system. However I know it is 'EPSG:4326'
I am writing a python script to perform tasks automatically, I don't want user interaction.

I am using the following commands :

fileName = "/path/to/raster/file.tif"
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(fileName, baseName)`

But the QGIS windows asked me the Coordinate reference system of this layer.
I tried to specified the crs while loading with something like

rlayer = QgsRasterLayer(filename=fileName, basename=baseName, crs='EPSG:4326')

But QGIS crashes.

Is it possible to specify a qml file where the crs is defined ?


texts: Changing a no value data from 0 to another value with python

texts:
I am using satellite data and by default it has a nodata value of 0 for each band.

This is fine except for when I'm making NDVIs where 0 is a real value we need to include and I don't want it to be classified as nodata.

I want to make sure all the bands/tifs I am working with have the same nodata value. So how do I convert a nodata value of 0 to something else, like -10000.

I was thinking of using arcpy's reclassify tool on my bands changing 0 to -10000 (before I make NDVIs) and then copying the raster with the copyraster tool using the setting for nodata as -10000.

This seems inefficient. Is there a tool that does this better? I'm open to open source libraries or arcpy


texts: How to add radius to map of over 32,000 zip codes

texts:
I'm incredibly new to creating these types of maps, but I'm trying to create a coverage map not unlike the cell companies show their coverage. Our database simply has lists of ZIP codes that our service covers and I want to display these on a map, colored in, like the cell phone companies do.

I have figured out how to do this on my own (using QGIS and joining my CSV list with the shapefile from the USPS of zip codes)... But the thing I can't figure out how to do is to show a radius around each zip coded area. Our coverage extends 100 miles (50 miles in each direction) from the center point of each zip code.

With a list of over 32,000 zip codes, is there a way to easily extend the coloring area around each zip code to include this extra area?


texts: Defining regions with certain size using polygon layer as input in ArcGIS for Desktop?

texts:
I'm working with ArcGIS 10.3 and have a project where I'm trying to define how may breeding habitats of a species that can fit into a specified area. I have a polygon layer with suitable habitat in this area and now I want to cluster all the suitable habitat polygons in this layer so it produces regions with 40 hectare suitable habitat (most polygons are smaller than that by themselves) but they do not need to be contacted to each other they can be spread out in an area up to 200 hectares.

I would like this to be optimized so as many territories as possible fit into this area. I'm thinking the region group tool can be used somehow if I convert the suitable habitat layer into raster?? Other methods are welcomed as well!


texts: Split tool is deleting segment at split point

texts: I turned editing on and selected this line feature from a feature class. After I selected this line feature, I clicked the Split Tool to separate the line at the vertex (the corner). When I click at the vertex, or any other point on the line, a portion of the line is simply disappearing instead of becoming two separate features. Any ideas on why this is happening?


texts:
I am using OSM basemap in my openlayers application.

    var OSMlayer = new OpenLayers.Layer.OSM( "Simple OSM Map");

    //Creating a Map
    map = new OpenLayers.Map('map', {
        maxExtent : OSMlayer.maxExtent,
        units : OSMlayer.units,
        numZoomLevels : OSMlayer.numZoomLevels,
        projection: new OpenLayers.Projection("EPSG:900913"),
        displayProjection: new OpenLayers.Projection("EPSG:900913"),
    });

    map.addLayers([OSMlayer]);

And I have a WMS service that using Web Mercator coordinates like (12.15, 88.65).

I am adding this WMS service on my map,

    layer = new OpenLayers.Layer.WMS("EcoRegions-Esri WMS", "http://myservicedomain.com/Regions/MapServer/WMSServer", {
            layers : "0",
            format : "image/png32",
            transparent : "true"
        }, {
            isBaseLayer : false,
            wrapDateLine : false
        });

    map.addLayers([layer]);

But my WMS tiles are not coming on my map.

when I look at chrome network, I see the WMS URL bbox are like this:

   myservicedomain.com/Regions/MapServer/WMSServer?LAYERS=0&
FORMAT=image%2Fpng32&TRANSPARENT=true&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&
STYLES=&SRS=EPSG%3A900913&BBOX=15028131.255,-5009377.085,20037508.34,0&
WIDTH=256&HEIGHT=256

But BBOX parameter should be like this:

BBOX=15.255,-50.085,20.34,0


texts: How to use st_geometry_type() to filter an sf object in r?

texts:
Put together this code to filter an sf object by geometry type as suggested by @Spacedman here

Posting here for others and for myself in case there are better ways to do this.

types <- st_geometry_type(my_sf_object)

types_df <- data.frame(types)

my_labelled_sf_object <- merge(my_sf_object, types_df,by.x=0, by.y=0, all.x=TRUE)

feature_type <-  "LINESTRING"

# value here can be any of
# GEOMETRY,POINT, LINESTRING, POLYGON, MULTIPOINT, MULTILINESTRING, MULTIPOLYGON
# GEOMETRYCOLLECTION, CIRCULARSTRING, COMPOUNDCURVE, CURVEPOLYGON, MULTICURVE,
# MULTISURFACE, CURVE, SURFACE, POLYHEDRALSURFACE, TIN, TRIANGLE

my_filtered_sf_object <- labelled_layers[labelled_layers$types %in% feature_type,]

this use of %in% seems correct but strikes me as very weird syntactically since really I'm selecting only rows with feature_type in that column rather than the other way around.


texts: Birch algorithm does not cluster as expected

texts:
I'm using Birch algorithm from scipy-learn Python package for clustering a set of points in one small city in sets of 10.

I use following code:

no = len(list_of_points)/10
brc = Birch(branching_factor=50, n_clusters=no, threshold=0.05,compute_labels=True)

In my idea, I would always end up with sets of 10 points. In my case now, I have 650 points to cluster, and n_clusters is 65.

But, my problem is that with too low threshold I end up with 1 address a cluster, just a tiny bigger threshold - 40 addresses per cluster.

What am I doing wrong here?


texts:
I have a global DEM and try to calculate contour polygons for 4 different heights in QGIS.
In the first image below, you can see the colored DEM (white = lowest; yellow = higher; orange = highest) as well as the overlying contour lines. The white area inside the yellow area is a sink.

After this i did convert the contour lines into polygons with the "Lines to polygons" function. Unfortunately, the sinks are not cut out in the polygon and remain the same class like the surrounding as you can see in the image below. The white areas from above have now the same class as the yellow area (its all pink).

Can someone explain me how I can archive, that the white area would have another class after polygonizing?


texts: QGIS 3 : Qt Signal emmitted when selecting a layer?

texts:
I am currently writing a QGIS 3 plugin. I would like to catch the signal emmited each time the user selects a new layer in the layer panel.

I have already searched in the QGIS API (QgsProject, QgsLayerTree, QgsLayerTreeView, ...)

QgsLayerTreeView has a signal currentLayerChanged but I cannot find how to connect QgsLayerTreeView to QgsProject.instance().

Does anyone has an idea on what is the solution or at least where I can find some information?


texts:
My first post here.
I'm designing a map at 500.000 for Italy and with EPSG 3035 or 3034 (Lambert Conical Conformal) I'm still deciding... anyway...
In my map I've several points whose symbol is a sort of "rose of winds" and I need they will remain aligned with the Grid North.
How could I do that? I think I'd have to calculate the angle between the Grid North and True North so I could put its value in a column that Qgis can read as an angle of rotation.

I tried with the formula

                        γ = arctan [tan (λ - λ0) × sin φ]

but I think there's something wrong coause the symbol doesn't align with the grid.


texts: Ghost Multipolygons block Overlay Intersection

texts:
I'm really desapointed, because I did everything to avoid Multipolygon in my script but last treatment blocks because of Multipolygon :

File "/.../Zonage/Zonage_standalone.py", line 254, in zonage
    zonage = gpd.overlay(clip, parc, how='intersection') # Intersection parcellaire et zonage

AttributeError: 'MultiPolygon' object has no attribute 'exterior'

You can find my script and datas here if you want to help me.


texts: How to issue http-request as soon a feature is clicked in ArcGIS map-viewer

texts:
I'm investigating ArcGIS 10.5 web client (map viewer, not sure what to call it). I've seen information how to populate a popup menu with data merged with template, or how to add custom actions to the popup by merging data from the feature.

We're looking for a way to issue an HTTP request as soon as a feature is clicked? The request being an URL (provided in the feature details) appended with the feature-id.

But we're not intending to customize the ArcGis map viewer application on a company level, because we're just one party out of many.

So my question is; is there built-in support to customize popup actions whenever specifically our layer is displayed? Or is the customization on a per viewer-app basis.


texts: GeoDjango, distance filter based on database column

texts:
I have a geodjango model with fields

coverage = 10 * 1609.344
reference_point = Point(-72,45) #or some other lat lng
class ServiceArea (geo_models.Model):
    center = geo_models.PointField(default=None, blank=True, null=True)
    radius = models.FloatField(default=floor(6 * 1609.344), null=True, blank=True)  # in meters

I am user center and radius to store a circle representing a service area.   I would like to query this to return all objects where the "center" is less than  (coverage + radius) away from "reference_point".

Or to put it another way, filter based on any part of the circlular service area being within a certain distance of a reference point.

This is simple with polygonfields, where I can do something like

ServiceArea.objects.filter(polygon__distance_lte=(reference_point, D(m=coverage)))

Although I haven't figured it out yet I imagine it would be possible to store the circles as polygons, and do the calculation as above, but this may have its own problems.  And it is simpler for other reasons to just store coordinates and radius in the database to approximate circles.

I've tried using django.db.models.F, but it doesn't seem to work in this situation.

What do you guys think?


texts: Reading table from geodatabase via ArcPy and doing what each row in field says after geopandas conversion?

texts:
Here is the table:

Reproducible example

import pandas as pd

d = {'id': [1, 2, 3], 'table_na': ['ER_ASTEN_ASTOT', 'INT_ASTEN_ASTOT','ER_ASTEN_ASTTOM']}
df = pd.DataFrame(data=d)
df

which gives:

   id   table_na
0   1   ER_ASTENOT_ASTOTA
1   2   INT_ASTENOT_ASTOTA
2   3   ER_ASTENOT_ASTTOM

I would like to read a table that looks like the above that I have in a gdb using Arcpy.

For every row in the table_na field, when it sees the ER or INT it should get the Erase or Intersection function in Arcpy respectively and use the relevant shapefiles (ASTENOT,ASTOTA,etc) for the calculation in a loop sequence.

I have already done this using geopandas but the process was very slow. In similar situations using ArcPy was very efficient. I do not have much experience in ArcPy which is why I am asking this question.

Here is the code in Geopandas I used to do the process and to help you help me:

#this reads the guide ex.ER_ASTENOT_ASTOTA from a txt.

import os
import fiona
import geopandas as gpd
rootdir = r'C:\Users\user\Desktop\a' # path to root directory you walk
sfiles = {} # a dictionary with all the .shp files
for entry in os.listdir(rootdir):
    dirpath = os.path.join(rootdir, entry)
    if os.path.isdir(dirpath):
        for file in os.listdir(dirpath): # Get all files in the subdirectories
            if file.endswith('.shp'): # If it's an .shp.
                filepath = os.path.join(dirpath, file)
                sfiles[filepath[:-4]] = gpd.read_file(filepath)

action_dict = {'INT': 'intersection', 'ER': 'difference'}
directory=input('Insert dir of the main folder')
with open(input()) as txtfile: #insert directory of txt
    x = txtfile.readlines()
for line in x:
    action, shape1, shape2 = line.split('_')
    print( action, shape1, shape2) # here line is ER_ASTOM_ASTIK or whatever line in your txt file
    if shape1 in sfiles and shape2 in sfiles:
        a=gpd.overlay(sfiles[shape1], sfiles[shape2], how=action_dict[action])


texts: I'd like to be able to pull up a DEM (or create a mosaic of them) and simply query the elevation at points I click. Is this possible using QGIS?


texts: Raster values in ArcMap need to be fixed for prediction maps

texts:
I am working on ArcMap with RGB rasters. The rasters have high values and i do not know how to fix it.
When i import them into R, it makes my prediction maps off. I have resampled them to another raster so they are in the same coordinate system and same pixel size.


texts: What prevents reducers from returning values?

texts:
I am attempting to calculate tree canopy coverage statistics (count, min, max, mean, median) on a census block level using the GLCF: Landsat Tree Cover Continuous Fields data set. I ran the script and discovered for some census blocks it will return all statistics, for others it will return only the mean and median, while yet other census blocks it won't return any of them. The code below references three census blocks that illustrate this behavior.

What causes this behavior?
Are there any changes I can make to guarantee all statistics are calculated?

_

var censusBlockIds = ['050534701002010',    // Returns all statistics
                      '051250105031020',    // Returns only mean and median statistics
                      '051250106006047'];   // Returns no statistics

// Get census blocks for analysis
var censusBlocks = ee.FeatureCollection("TIGER/2010/Blocks")
                     .filter(ee.Filter.inList('blockid10', censusBlockIds));

// Get 2010 Landsat Tree Cover Continuous Fields
var treeCanopyCover = ee.ImageCollection('GLCF/GLS_TCC')
                        .filter(ee.Filter.date('2010-01-01', '2010-12-31'))
                        .filterBounds(censusBlocks)
                        .select('tree_canopy_cover')
                        .max();

// Function to calculates statistics and add block id
function censusBlockTreeCoverage(censusBlock) {
  return  calcTreeCoverage(censusBlock.geometry())
    .set('blockid10', censusBlock.get('blockid10'));
}

// Function to calculate tree canopy coverage statistics for census block
function calcTreeCoverage(censusBlockGeometry) {

  // Count pixels census block intersects
  var countDict = treeCanopyCover.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: censusBlockGeometry,
    scale: 30,
    maxPixels: 1e9,
    bestEffort: true
  });

  // Calculate minimum tree canopy cover for census block
  var minDict = treeCanopyCover.reduceRegion({
    reducer: ee.Reducer.min(),
    geometry: censusBlockGeometry,
    scale: 30,
    maxPixels: 1e9,
    bestEffort: true
  });

  // Calculate maximum tree canopy cover for census block
  var maxDict = treeCanopyCover.reduceRegion({
    reducer: ee.Reducer.max(),
    geometry: censusBlockGeometry,
    scale: 30,
    maxPixels: 1e9,
    bestEffort: true
  });

  // Calculate mean tree canopy cover for census block
  var meanDict = treeCanopyCover.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: censusBlockGeometry,
    scale: 30,
    maxPixels: 1e9,
    bestEffort: true
  });

  // Calculate median tree canopy cover for census block
  var medianDict = treeCanopyCover.reduceRegion({
    reducer: ee.Reducer.median(),
    geometry: censusBlockGeometry,
    scale: 30,
    maxPixels: 1e9,
    bestEffort: true
  });

  // Return new geometry-less feature with calculations
  return ee.Feature(null)
    .set('treeCoverageCount', countDict.get('tree_canopy_cover'))
    .set('treeCoverageMin', minDict.get('tree_canopy_cover'))
    .set('treeCoverageMax', maxDict.get('tree_canopy_cover'))
    .set('treeCoverageMean', meanDict.get('tree_canopy_cover'))
    .set('treeCoverageMedian', medianDict.get('tree_canopy_cover'));
}

// Define visualization parameters
var treeCanopyCoverVis = {
  min: 0.0,
  max: 100.0,
  palette: ['ffffff', 'afce56', '5f9c00', '0e6a00', '003800'],
};

// Calculate statistics for each census block
var censusBlockStatistics = censusBlocks.map(censusBlockTreeCoverage);
print(censusBlockStatistics);

// Add layers to map for visual confirmation
Map.addLayer(treeCanopyCover, treeCanopyCoverVis, 'Tree Canopy Cover');
Map.addLayer(censusBlocks, {}, 'Census Blocks');
Map.centerObject(censusBlocks);


texts: import arcpy yields "TypeError: This object does not support enumeration"

texts:
Similar to Enumeration Problem with Eclipse and arcgisscripting?, with Arcgis v10 I no longer seem to be able to run python scripts from a command shell. Entering the same code into an interactive python console from the same shell does work! Other scripts I have from a few months ago seem to work fine. I don't get it.

This works:

d:\s>python
Python 2.6.6 (r266:84297, Aug 24 2010, 18:46:32) [MSC v.1500 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import arcpy
>>> arcpy.env.workspace = 'd:/'
>>> print arcpy.env.workspace
d:/
>>>

Create a minimal test file:

d:\s>copy con d:\code\arcgis\test-arcpy.py
import arcpy
arcpy.env.workspace = 'd:/'
print "Workspace is: ", arcpy.env.workspace^Z
        1 file(s) copied.

and give it a spin:

d:\s>python d:\code\arcgis\test-arcpy.py
Traceback (most recent call last):
  File "d:\code\arcgis\test-arcpy.py", line 1, in <module>
    import arcpy
  File "C:\ESRI\ArcGIS\Desktop10.0\arcpy\arcpy\__init__.py", line 17, in <module>
    from geoprocessing import gp
  File "C:\ESRI\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing\__init__.py", line 14, in <module>
    from _base import *
  File "C:\ESRI\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing\_base.py", line 568, in <module>
    env = GPEnvironments(gp)
  File "C:\ESRI\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing\_base.py", line 565, in GPEnvironments
    return GPEnvironment(geoprocessor)
  File "C:\ESRI\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing\_base.py", line 521, in __init__
    self._refresh()
  File "C:\ESRI\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing\_base.py", line 523, in _refresh
    envset = (set(env for env in self._gp.listEnvironments()))
  File "C:\ESRI\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing\_base.py", line 523, in <genexpr>
    envset = (set(env for env in self._gp.listEnvironments()))
  File "C:\Python26\lib\site-packages\win32com\client\dynamic.py", line 246, in __getitem__
    raise TypeError("This object does not support enumeration")
TypeError: This object does not support enumeration

update it's directory specific! If I copy that same test-arcpy.py file somewhere else it works.


texts: Public Geoserver preview via third-party website?

texts: A few years ago I remember seeing a website that would connect to a public Geoserver offering and allow you to preview all available layers overlaid on a Google Maps or similar base layer.  Does anyone else remember that?  Is it still around somewhere?


texts:
Is there is a way in QGIS (or SQL, ArcMap) to split multi-geometries in polygons automatically. I could do this manually (split features, split multi-part to single-part etc...), but this would just be so time consuming. Here's a visual of what I'm wanting the polygons to look like:

As you can see polygon 2 has been broken up into 5 polygons instead of being a single polygon in polygon 1. I tried what this person has done in SQL HERE, but then that would only break polygon 1 into three polygons still having a multi-geometry polygon.

Is there a way to automate this process?


texts: Print composer: Frozen layers disappeared after server migration

texts:
I am currently working on a project file with around 30 print composers using different layers. Some of the layers come from a PostGIS database connection, others were imported through the Add Layer menu. Some layers have styles applied from locally saved .qml files, other styles were set manually. Every print composer has at least one map window and every single one of those map windows has its layers frozen (so I can modify the styles and the layer order for one of the composers without affecting the 30 other composers).

Now, I work in an office setting and the project, as well as some layers and all of the .qml files, are saved on a network mapped drive. The address of the mapped drive was recently changed due to a migration on a newer server. The first time I opened one of the print composers after the server migration, I noticed that while all the layers did appear in the project itself, none of them would show up in the print composer as they should, considering that were frozen in a map window. I assume that the Freeze Layers function depends on an absolute path to the layers or something.

This is a major inconvenience as it means that I will have to manually re-arrange the layer order, styles and labeling for 30+ print composers. Does anyone have an explanation as to why this happened? Is there a way to fix it, like by fiddling with the .qgs file?


texts: Given list of points (lat, long), how to find all points within radius of a given point

texts:
It's my first time working with GIS. I really wish I could share some code, but I don't even know how and where to begin. I have a Pandas dataframe with around 10,000 records. Each of them have a latitude and longitude.
Assume that after applying some filters, I end up with 10 records relevant to me.

Now, I want to find all the records which have coordinates lying within 20 miles (radius...?) of ANY of these 10 records.

The dumbest brute force solution I thought of was using geopy and looping over the records and finding the distance. Then I came across this solution to use kdtree - https://stackoverflow.com/a/41155458. But the comment below it says that it isn't accurate.

And I also tried to look at geopandas, and use the buffer and unary_union functions - http://geopandas.org/geometric_manipulations.html. The only thing I have with me right now is the dataframe with the latitudes and longitudes.

How should I proceed?


texts: Map service appears in the web REST view but not in ArcMap?

texts:
I have an ArcSDE server serving some poly layers through a map service running on ArcGIS Server 10.0.  The service appears to be working--I can navigate to it just fine in the REST view and I can see the data when I do the View In ArcGIS JavaScript option.

However when I load up ArcMap and attempt to add the service to a layer it appears to succeed but does not show up; the screen remains white no matter what extent I'm viewing.  The layer is set to be viewable at all extents and appears to be the correct spatial reference so I'm a bit confused why the javascript viewer can see it but not ArcMap.

edit I should add that I have a Silverlight app pointing at the service as well and that also does not show the data.


texts: Is it possible to set a layer style (e.g. set the color ramp) before loading the layer using pyQGIS?

texts:
I am implementing a tool which allows an user to calculate some statistics with some data from the registration office. My question is, can I set some layer styles using pyQGIS?

The analysis part is already finished and in a next step I want to laod the layer (which also works perfectly) and visualise it in QGIS. The problem is that I have found until now no way to "say" QGIS that the layer should be visualised by using "graduated colors", a special color ramp and a special classification.

Is there a way to do this using pyQGIS (mayber by using map canavas) or does the user has to do it manually in the layer settings?

Regards from Hamburg,


texts: Adjusting Polygon Visibility with CartoCSS

texts:
I am attempting to adjust the visibility of a polygon layer in Carto by zoom level.  To keep it simple for now I am looking to have the opacity of 1 from zoom 0-7 and opacity 0 from zoom 8+.  I am working directly in Carto (not tilemill, or elsewhere) and have tried a few methods.  Most notably, this method (found on this question) only worked for point data, not polygons.

Getting the zooms to adjust for point data worked using the following, now I'm just trying to adapt it for use with polygons.

    #layer {
    marker-allow-overlap: true;
    [zoom = 0]{opacity:1}
    [zoom = 8]{opacity:0}
    }

Right now I've gotten myself stuck at

    #layer {
    [zoom = 0]{opacity:1}
    [zoom = 8]{opacity:0}
    }


texts: ATAK connection to OS WMTS raster server requests username and password

texts: Attempting to connect to the OS raster map WMTS server from WinTAK asks for a username and password. Trying my OS account email address and password fails. Connection using same link works fine in QGIS. Has anyone encountered anything similar before?


texts: OpenLayers is not defined while accessing over the html and js

texts:
The above is my working directory, all html files are under WebContent folder and the problem is that when all the code is written in a single html file working fine.
If I split the OpenLayers code into 'js' files under 'js' folder. The files within 'js' folder couldn't refer OpenLayers library(OpenLayers.js).

ReferenceError: OpenLayers is not defined
So that var's not getting initialized.

{

        <meta name="viewport" content="width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=0">
        <link rel="stylesheet" type="text/css" href="ext-321/resources/css/ext-all.css" />

        <link rel="stylesheet" type="text/css" href="css/style.css" />

}

working with localhost.
and calling init() method explicitly from Ext.onLoad() while I am calling with in the same page it works fine and there is no problem with map it is displaying

Ext.onReady(function(){
        Ext.QuickTips.init();
        Ext.BLANK_IMAGE_URL = "ext-321/resources/images/default/s.gif";
        // turn on validation errors beside the field globally
        Ext.form.Field.prototype.msgTarget = 'side';
        new Ext.Viewport({
        layout      :   'border',
        renderTo    :   Ext.getBody(),
        region  :   "center",
        //margins   :   '0 5 5 0',
        defaults:   {border:false},
        items       :   [{
            region  :   "center",
            layout  :   "fit",
            id      :   "center_region",
            /*//tbar    :   new Ext.Toolbar({//551
                            id: 'map-toolbar'
                        }),
            html    :   "

"
            *///555//
            items   :[{
                xtype   :"tabpanel",
                id      :"tab_map_status",
                activeTab: 0,
                defaults    :   {border:false},
                items:[{
                            title   :   'Map',
                            id      :   'tab-map',
                            tbar    :   new Ext.Toolbar({
                                            id: 'map-toolbar',
                                            items: []
                                        }),
                            html    :   "

"
                        }]}]});init();});'

Please help me.


texts: How does QGIS Topology Checker check the topology of SHP and OSM files?

texts:
Can anybody make me clear if .shp & .osm (OpenStreetMap Data) are holding topological data/relations or not.
if .shp files are not holding then how qgis plugin topology checker checks for the topology of these files.


texts: Choosing pixel size when projecting raster from geographic to projected coordinate system

texts:
I have a raster showing animal abundance (straight count) in all of Canada. The raster was created unprojected in GRS 80 with a pixel size of 0.0833 decimal degrees. I was ask to project this raster to Canada Albers Equal Area Conic. I have never seen an unprojected raster before and don't know why it would be created this way. When inputting the raster into ArcGIS 10.3 reproject tool, the default pixel size is set to 8787.96 m. I realize that no one pixel size would be correct because a decimal degree in one part of the country would be a different length than another part.

Is there any standard protocol or guidelines in this situation?


texts:
I am trying to share a map as a service and overwrite an existing Feature Service in ArcGIS Online. I put together a small model to connect to the database, to publish and replace existing item online, but I am not seeing the "Share As a Service" tool in the ArcToolbox.
How can I share it as a service with ModelBuilder?


texts: Scripting (Python) Feature Class content into Shapefiles following an UTM data structure

texts:
Consider a Feature Class and two of its fields are "UTM_grid" and the "OID".

Scripting with python 2.7 for ArcGIS 10.3.

"UTM_grid" field contains the geographic info (e.g. B1SD40B) and was broken down into "UTM_block" + "UTM_sheet" (B1 + SD40B) and then using a dictionary and list structure UTM_block-->UTM_sheet-->OID was organized with this syntax:

{ 'UTM_block_1' : { 'UTM_sheet_A' : [ oid_1 , oid_2 ] }, { UTM_sheet_2 : [ oid_3 , oid_4, oid_5 ] } }

We need to create shapefiles with features (identified by the object value - oid) based on the follow data structure (OIDs will be organized regarding an UTM/location based structure):

UTM_block_* will be a folder
UTM-sheet_* its subfolder.
List of objects (OIDs) will populate the shapefile.

Note:
Features Classes to export are in a SDE geodatabase. We read/list the FCs in that GDB and for each FC we need to process the export of its features by UTM block-sheet into shapefiles. Looking at the fields "UTM_grid" and "ObjectID" we check while sheet ( utm_grid[1] ) = "foo", read the corresponding OID, and put the Feature into a Shapefile inside the (sub)folder (which has the associated sheet name). The same logic could be used to organize sheets inside the block folders.

Another point we have to deal with is that the shapefile should only have some fields of the original Feature Class:

We have used fieldmappings and arcpy.FeatureClassToFeatureClass_conversion() to create a non-UTM-organized shapefile from a Feature Class with the fields that we need.

Now the problem is how to process the feature export to get UTM-organized shapefiles! How can we manage it?

Bellow is what we have now inside a function:

with arcpy.da.SearchCursor(feature_class, ['UTM_grid' , 'OBJECTID']) as cursor:
    dic = dict()
    for row in cursor:
        UTM_sheet_value = row[0]
        oid_value = row[1]

        try:
            b_s_split = UTM_sheet_value.split('_') # split of block & sheet

        except Exception, e:
            pass
            print(UTM_sheet_value)
            print(e)

        else:
            dic.setdefault( b_s_split[0], {} ).setdefault( b_s_split[1], [] ).append( oid_value )

            # create folders & subfolders (blocks & sheets) in a directory
            dir_path = os.path.dirname("E:\\")
            dirs = [ [b_s_split[0]] , [b_s_split[1]] ]

            for item in itertools.product(*dirs):
                if not os.path.isdir(os.path.join(dir_path, *item)):
                    os.makedirs(os.path.join(dir_path, *item))


texts: What is the Python expression for calculating Lat/Long in Decimal Degrees in ArcGIS?

texts:
I'm calculating 4 double fields for a polygon feature class (exploded, single-part features): SquareFeet, Acreage, Lat, and Long

I found the correct python expressions to calculate SquareFeet & Acreage using the Calculate Field (Data Management) tool.

For the Lat and Long fields, I've tried using !SHAPE.CENTROID.Y! which calculates in either feet or meters depending on which unit is set on the data, but since I need the calculation in decimal degrees this isn't acceptable. And I've tried !SHAPE.CENTROID.Y@DECIMALDEGREES! to get the calculation in decimal degrees, but this expression just flat out doesn't return anything. I've tried setting the parser to PYTHON and PYTHON_9.3.

Does anyone know if there's a python expression to calculate latitude and longitude in decimal degrees?

(!SHAPE.CENTROID.X@DECIMALDEGREES! of course was used to try and calculate Long)


texts: Drawing line with specified length using QGIS?

texts: How to draw a line with specified length using QGIS 2.16.1?


texts: Compare planet downloaded images with different resolution in Python

texts:
I have downloaded satellite images using Planet API for an Area of Interest. (Image Details: PSScene4Band, analytic)The problem is Planet returns images with different number of pixels for the same subarea, with varying dates.
I want to compare NDVI pixel by pixel of these images. Since the number of pixels differs, it's not possible.

Attached images I downloaded.Images link here

File 20180502_181704_1025_subarea.tif has size: 2027, 1193.

File 20170919_180930_0f22_subarea.tif has size: 2028, 1194

My code

import numpy, sys
from osgeo import gdal
from osgeo.gdalconst import *
import numpy as np
import random
import matplotlib.pyplot as plt

# register all of the GDAL drivers
gdal.AllRegister()

# open the image 1
inDs = gdal.Open("20170919_180930_0f22_subarea.tif")
if inDs is None:
   print('Could not open image file')
   sys.exit(1)
r=np.array(inDs.GetRasterBand(3).ReadAsArray(), dtype=float)
n=np.array(inDs.GetRasterBand(4).ReadAsArray(), dtype=float)

np.seterr(divide='ignore', invalid='ignore')
#Ignore the divided by zero or
Nan appears
ndvi1=(n-r)/(n+r) # The NDVI formula

# open the image 2
inDs2 = gdal.Open("sample/20180502_181704_1025_subarea.tif")
if inDs2 is None:
   print('Could not open image file')
   sys.exit(1)

r=np.array(inDs2.GetRasterBand(3).ReadAsArray(), dtype=float)
n=np.array(inDs2.GetRasterBand(4).ReadAsArray(), dtype=float)

np.seterr(divide='ignore', invalid='ignore') #Ignore the divided by zero or
Nan appears
ndvi2=(n-r)/(n+r) # The NDVI formula

ndvi = ndvi2 - ndvi1

Problem occurs in the line ndvi = ndvi2 - ndvi1
 since the numpy array dimensions differ, as the image resolution (rows x columns) differs.
Note: I assume ignoring the pixels may result in the wrong result; as I have no clue which pixel is missing or extra.

I did make sure that both the TIFF images have same corner coordinates, but pixels differ.


texts: Assigning elevation to point shapefile based on point file using ArcGIS Desktop

texts:
Goal:

Extract elevation data (Data B) based on given point file locations (Data A).

Current problem:

I have several sets of xyz data for large areas with a resolution of 1m. I converted those to Point Data (Data B). When I apply the approach of (1) Point to Raster, (2) Extract Points from Raster it requires a lot of time for the geoprocessing. That's why I'm searching for a way to bypass the rasterprocessing.

I was wondering if its possible to just assign the elevation of the closest points (Data B) to my target points (Data A)?

I found some similar questions, but none of them discussing which method is the fastest and most of them were handled in QGIS, while I work with ArcGIS 10.7:

Assigning a nearest polygon to a point
Joining attributes from nearest point in QGIS?
Assigning Elevation data from one point shapefile to another point shapefile in QGIS


texts:
I followed the instructions for loading OSM data into pgrouting found here: http://anitagraser.com/2011/12/18/osm2po-part-2-pgrouting-on-osm-the-easy-way/

Everything works fine using osm2po to calculate a route, so I assume I setup everything correctly. However...

In QGIS when I load the database and fire up the pgrouting layer plugin, I get an error when I try to select the node nearest node to calculate a route from.

"Relations Roads Does not exist:

Any ideas or help?


texts: How to create Geodesic Polygons in OpenLayers?

texts:
I would like to create geodesic polygons in openlayers. This is not the same as this question although VERY similar.

Geographika was able to do it but his solution is based on a circle. I'm interested in creating polygons that aren't necessarily circles, such as squares. Or maybe even a line.

Anyone know of any javascript snippet that might convert my geometric openlayers polygons to geodesic ones?

Google Maps has the Geodesic:true property for its shapes. OpenLayers should definitely have it.


texts:
I am trying to run a multiprocessing on a zonal statistic in pure Python. I am new on it.  I have to apply a zonal statistic with a shapefile with more than 6 million features on a taster of 6 GBs (dimension X: 639760 Y: 452420 ).
I am implemented a multiprocessing starting from this script this script (from @Perrygeo!!) but the process run forever; after two weeks it haven't finished yet.

The process was reading the entire shapefile and the entire raster and after this processing using the code.

So I would like to implement this process:

reading the shapefile piece by piece from a PostGIS database. I am not able to read a shapefile piece by piece with geopandas or fiona.
applying the zonal statistic only on the part intersecting the raster, a very big raster (6.3 GB also with LZW compression). I am assuming the processing will be done only on the intersection area.
printing the final output on the database

Can I use asyncio to implement a concurrency to speed the process and applying at the same time a multiprocessing on the zonal statistic itself?
I need to write a fast code...


texts: How to rasterize point data into a grid which can ensure 1 point in each cell at most

texts:
I'm trying to create a raster whose grid cell takes at most one point. The data looks like this:

> head(ras)
       lon    lat         z1    z2
1 267573.9 2633781 213.29545     6
2 262224.4 2643701  69.78261    15
3 263742.7 2670841  51.21951     1
4 259328.4 2739781 301.98413    29
5 264109.8 2463763 141.72414    11
6 255094.8 2063428  88.90244    35

z1/z2 are measurements of two variables, which can be put into separate layers, i.e. layer 1 takes z1 values, and layer 2 takes z2 values. The two layers use the the same grid system. Empty cells can be assigned with 0 or NA.


texts: Any Open Street Map SDK or library for Windows CE 6.0

texts: I have been assigned a task to develop a navigation app for Windows CE 6.0 and I was going through Frameworks and there is only osmsharp that supports windows phones, are there any sdk or libraries specifically for Windows CE 6 development that fullfill the need of development of a navigation app at least having basic features along with offline?


texts: QGIS set all layers with custom transparency with python?

texts:
I'm exploring QGIS python functions to set all raster layers with a custom transparency, i.e 0,0,0,100 .

Looking at some previous questions. There is a way to loop all the layers, but the aspect to set transparency value on RGB values, I've been having difficulties understanding the documentation.

How to set the transparency for multiple layers or add a global transparancy preference?

for layer in iface.legendInterface().layers():
   layer.renderer().setOpacity(0.2)

How do I set layer transparency in QGIS 2.0 with Python?
The example here uses TransparentSingleValuePixel but i notice there is a TransparentThreevaluePixel

x = QgsRasterTransparency.TransparentSingleValuePixel()
x.pixelValue = 123
x.transparencyPercent = 50
rt.setTransparentSingleValuePixelList( [ x ] )


texts: Importing multiple Esri ASC files into GRASS, then exporting to Unreal Engine

texts:
I'm trying to import LiDAR scan data from Defra in the UK into Grass.  The data appears to be in Esri ASC format, in multiple separate files.

As this is the first time I've used GRASS I'm struggling a bit!  I've managed to get one tile in using the Grass GUI and Import raster data - Import of common raster formats r.in.gdal.  The issue is I can't work out how to import multiple tiles into the same map.

The reason I'm trying to do this is to then create something I can import into Unreal Engine 4.


texts:
Downloading data via the new soilgrids.org is working nicely for all layered information, but not for WRB classification. I tested it at different PC and different browsers, but it always gives this error message:

<ows:ExceptionReport version="2.0.1" xml:lang="en-US" xsi:schemaLocation="http://www.opengis.net/ows/2.0 http://schemas.opengis.net/ows/2.0/owsExceptionReport.xsd">
<ows:Exception exceptionCode="ExtentError" locator="extent">
<ows:ExceptionText>
msWCSGetCoverage20(): WCS server error. Raster size out of range, width and height of resulting coverage must be no more than MAXSIZE=8192
</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

The request being made is:

https://maps.isric.org/mapserv?map=/map/wrb.map&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=MostProbable&FORMAT=image/tiff&SUBSET=long(0,2)&SUBSET=lat(0,2)&SUBSETTINGCRS=http://www.opengis.net/def/crs/EPSG/0/4326&OUTPUTCRS=http://www.opengis.net/def/crs/EPSG/0/152160

Would be interesting to know if this is an error due to my hard-/software or a general problem (and how to solve it).


texts:
I have polygon feature and want to be able to generate points inside it. I need this for one classification task.

Generating random points until one is inside the polygon wouldn't work because it's really unpredictable the time it takes.


texts: Exporting layer to shapefile after styling using Python in QGIS

texts:
I'm trying to use Python to add a GML to my QGIS 3.4 project, apply styling to it, and then export the layer as a shapefile.

It exports the layer, but it does so before applying the styling, which uses  Rule-based symbology, which filters some data away that I want separated in different shapefiles.

I think it has to do with layer it exports, but I need it to export the layer after it added the styling to it,

writer5 = QgsVectorFileWriter.writeAsVectorFormat(layer5,directory + shpmap + "/road1.shp","utf-8",driverName="ESRI Shapefile")

Does anyone have an idea how to do this?

import os

from qgis.core import QgsProject
directory = QgsProject.instance().readPath("./")
print(directory)
root = QgsProject.instance().layerTreeRoot()

styles = "/Styles"
importfolder = "/Import"
shpfolder = "/Shpexport"
dxffolder = "/Dxfexport"

group1 ="Group"
group2 = "Analyse"
group3 = "Analyse2"

file_to_open5 = directory + importfolder +"/road1.gml"
file_to_open6 = directory + importfolder + "/road2.gml"

style5 = directory + styles +"/road1.qml"
style6 = directory + styles + "/road2.qml"

group1 = root.addGroup(group1)
group2 = root.addGroup(group2)
group3 = root.addGroup(group3)

layer5 = QgsVectorLayer(file_to_open5, "road1", "ogr")
if not layer5.isValid():
    print ("road1 failed to load!")
else:
    print ("road1 was loaded successfully!")
    layer5.loadNamedStyle(style5)
group1.addLayer(layer5)

layer6 = QgsVectorLayer(file_to_open5, "road2", "ogr")
if not layer6.isValid():
    print ("road2 failed to load!")
else:
    print ("road2 was loaded successfully!")
    layer6.loadNamedStyle(style6)
group1.addLayer(layer6)

writer5 = QgsVectorFileWriter.writeAsVectorFormat(layer5,directory + shpmap + "/road1.shp","utf-8",driverName="ESRI Shapefile")
writer6 = QgsVectorFileWriter.writeAsVectorFormat(layer6,directory + shpmap + "/road2.shp","utf-8",driverName="ESRI Shapefile")


texts: How to find out to which road each points belongs to?

texts: I am new to GIS . I hava a folder it include a .mxd file and a .gdb folder, the .gdb folder hava many .gdbtable and .gdbindexes  files etc, it's a road network data of a city. Now, I want to know which road it belongs to for a given latitude and longitude, which software or sdk should I use? is there have a similar open source project?


texts: Hec-Ras 6.0 error creating terrain (Step 1 Translating to GeoTIFF)

texts:
I want to create a terrain in RAS Mapper (right click 'Terrains' -> 'Create a new RAS Terrain'). Following error occure:
Computation Task                                                                 |   hh:mm:ss
---------------------------------------------------------------------------------------------
Importing 1 of 1: dgm_sa_hec.asc -> Terrain.dgm_sa_hec.tif
Step 1 of 4: Translating to GeoTiff with SRS...
ERROR 1: missing [

ERROR 1: Failed to process SRS definition: ESRI::*path*\DGM\dgm_sa_hec.prj
Step 1 of 4: Translating to GeoTiff with SRS...                                  |          0
Error creating Terrain: *path*HecRas\Terrain\Terrain.dgm_sa_hec.tif: No such file or directory
Cleaning up files...

I already used an ASCII File and a TIFF but both doesn't work. There are no special characters in the path. The Terrain folder is created automaticly.
Does anyone have experience with a similar problem?
Using hec-ras 6.0 Beta2


texts: Alternatives to GeoRSS for time-sensitive point data

texts:
I'm looking for a suitable data format for various organisations to publish point data with a recent time component, like events, road closures, etc.

GeoRSS is the obvious candidate but:

It's XML, which is a bit painful, hard on the eyes etc
It's GML, which is not as widely supported as, say, GeoJSON
Because it supports more than just points, it's more verbose than needed for this case.

I'm tempted to just use a reverse-sorted CSV file:

date,id,lat,lon,description

Is there any reason a "feed" of data with timestamps needs to use a different format than any other geospatial data, anyway?


texts: Something wrong with The size of the converted shapefile (arcpy.PolygonToRaster_conversion)

texts:
I currently have a Python script to process some vegetation data, while using some specific shapefiles.

The first step is to build a geo lookup table which needs to convert the shapefile from polygon to raster and then from raster to ascii. And the extent of the converted file must match the raw data file.

Even though I set up the environment variables based on the raw data file, the result still show a different extent with the raw data file. Normally they should match since it did work when I tested it with ArcMap manually. However, it just did not work with ArcPy and Python code. Any opinions?

Below is the snippet of the code:

env.workspace = workingdir
env.overwriteOutput = True
arcpy.env.coincidentPoints="MAX"
arcpy.env.extent = rawdatafile
arcpy.env.cellSize = rawdatafile
arcpy.CheckOutExtension("Spatial")
arcpy.CheckOutExtension("GeoStats")

arcpy.PolygonToRaster_conversion(in_features=shapefile, value_field="FIPS", out_rasterdataset=shpraster, cell_assignment="CELL_CENTER",cellsize=rawdatafile)
arcpy.RasterToASCII_conversion(in_raster=shpraster, out_ascii_file=shpascii)`

arcpy.RasterToASCII_conversion(in_raster=rawdatafile,out_ascii_file=dataAscii)

However, the dataAscii and shpascii have different numbers of columns and rows.


texts:
I have a couple of thousand polygons in SpatiaLite.  I am trying to do a "touches" query:

select map1.* from map1,map2
where touches(map1."Geometry",map2."Geometry")

and wow, is it SLOW!

However, if I ask it to just do it for one parcel in map1, it runs really fast.

select map1.* from map1,map2
where touches(map1."Geometry",map2."Geometry")
and map1."ROWID" = 753

I do expect that the first query will run slower, but it is amazingly slow.  It runs very fast in SQLServer, Manifold GIS, and PostGIS.  Is Spatialite just really inefficient?


texts:
What are the best books/web sites for C# development with ArcObjects?

The ArcGIS Resource Center is very helpful, but I am trying to find sources with more examples.


texts: Reading lines within Bounding box failed

texts:
I am developing an application to display OSM data. In my scenario I need to display the data for one entire country. So my approach is that I will load curresponding countries osm data to postGreSQL and fetch it back on demand.

I am providing the zoom facility here (from 0 - 19) and for each zoom  level I will query different data based on what need to be displayed at that zoom level and using a bounding box. I am new to this approach. As there exist experts here, I would like to know how feasilible would be this approach in terms of performance? (Ealrlier I try to use osm.pbf files and trying to read it dynamically and display; but it hit the performance badly.

My second question is: I tried to query all the lines within a given bounding box but it displays empty row as result. May I know what I am doing here now?

SELECT * FROM public.planet_osm_line
where ST_INTERSECTS(way ,ST_MakeEnvelope(1.4411000, 41, 1.4411530, 45));

Neither the below query worked (I assume that I can give the lattitudes and longitudes here to define the bounding box here reference: Select bounding box using postgis)

SELECT * FROM public.planet_osm_line
where ST_INTERSECTS(way ,ST_MakeEnvelope(-180, -180, 180, 180));

I tried the below query also to fetch some data in France area. But failed again

SELECT * FROM public.planet_osm_line where ST_INTERSECTS(ST_Transform(way,4326) ,ST_MakeEnvelope(-73.3968207647942, 49.296681801433,-73.823777534737, 49.5751193739016));

Please advise.


texts: Use OS-specific file dialog to relink layer when file moves?

texts:
When I move a file out from under QGIS, of course I have to tell the program where it moved to.
If you close a QGIS project and then move or rename a layer source file, you get a prompt when you reopen the project asking where you moved the file, and you can use a OS-native (in my case, MacOS-native) open-file dialog to indicate where the source file moved to.
However, when I move a source file while QGIS is open, I can't figure out how to get an open-file dialog to link it to the changed source file. All I get is a clunky built-in file navigator which is not intuitive to use.
My default now is to just close the QGIS project and reopen it when I move layer source files, but is there a better way?


texts:
Using the topological corrector inside Qgis I have identified overlapping polygons. My situation is represented into this image

My polygon of interest is located below another polygon (the red polygon). I would like to correct this error, however I can't select the red polygon to delete.

Do you have any other idea to resolve this situation?


texts: Repeating Python script with output from previous step

texts:
I have a raster with gaps in it and want to fill the gaps with the mean of the surrounding cells using focal statistics. I found the code to achieve this here
Fill in NoData gaps in raster using ArcGIS for Desktop?

I then want to take the output from the focal statistics and repeat about 15 times in order to expand the raster.

I have managed to get it to work somewhat but it stops after about 3 iterations

arcpy.env.workspace= r'C:\Users\Environmental\Tides\Analysis\Analysis.gdb'

arcpy.env.overwriteOutput = True

Waves = arcpy.sa.Con(arcpy.sa.IsNull("M2_reproj.tif"),arcpy.sa.FocalStatistics("M2_reproj.tif",arcpy.sa.NbrCircle(1),'MEAN'), "M2_reproj.tif")
Waves.save("C:\Users\Environmental\Tides\Analysis\Waves.tif")

seq = [1,2,3,4,5,6,7,8,9,10,11,12,13,14]

for num in seq:
    Waves = arcpy.sa.Con(arcpy.sa.IsNull("Waves"),arcpy.sa.FocalStatistics("Waves",arcpy.sa.NbrCircle(1),'MEAN'), "Waves")
    Waves.save("C:\Users\jtwh3\TOM_GIS\Environmental\Tides\Analysis\Waves.tif")


texts:
I would like to hide/remove the checkbox of a LayerNode, since I need to have it always on display.

Is there any sort of option to do this, or must I use the Ext object hierarchy to locate the checkbox object and remove it?


texts: Joining Textfile to Shapefile Arc 10 and 10.2 Differences

texts:
I am trying to join a textfile generated in stata that includes caoname, caoid, and fips. The fips code should match to a shapefile county map that matches to the cntyidfp category.
In ArcGIS 10.0, doing a straightforward join based on these attributes works. However, in 10.2, ArcGIS only allows for the joining of the cntyidfp to the caoname, which is completely useless.
What's the difference between the two versions that makes this difficult?
I have tried exporting to an excel document first, ensuring both the fips and cntyidfp data have the same number of characters, and creating a new table of the textfile within arcgis, but none of these solutions seem to work.
What am I missing?


texts: QGIS: Join attributes by location calculates values wrong

texts:
It seems like I got the same problem which is discussed in this thread. But since it's 2 years old already, there's an easy solution to this by now?
I'm joining a vector layer which contains seperat buildings with values to a grid with the algorithm "join attributes by location (summary)". My "summaries to calculate" option is "sum".
Now my problem is, sometimes the buildings cross the line of a rectangle of the grid (see image). In these cases, the building gets recognized twice and thus the resulting sum of the layer is higher than it should be.
My question is, how can I tell QGIS to put the buildings each only in one grid if they're set in two?
The best solution would be, if the buildings would automatically be set in the square where the most of it is in.


texts: Batch distance calculation with v.net.distance

texts:
I have three shapefiles:

one containing destinations
the second containing postcodes
the third containing the road network

I want to calculate the distance between the postcodes and destinations with the help of GRASS v.net.distance in QGISS 3.2 GUI.

I can easily calculate the distance with this tool between two points only but not between multiple feature(red dots show the destinations, the pink ones show the post codes).

Can anyone explain how to run the batch process for this?


texts:
i like to get the neighbors of the Parcels "P1" like following :

From B1 to B2 => P5

B2 => P4

From B2 to B3 => P3

From B3 to B4 => P2

From B4 to B1 => P6

i like to make this process automatic by a GIS software like ArcGIS or Postgis, i have thought to make a buffer of the parcel P1 , to get the lines of the boundary of P1 but i will not get the parcel P4 as neighbor .

any solution to get the neighbors like described above.


texts:
i'm having an issue with QGIS 2.14.3 (Windows 10). I've two layers:

PostGIS table (EPSG:23030)
Shapefile (EPSG:25830)

I setup on the fly transformation (23030) in QGIS, just for check, but both layers does not overlay.

However if I do same thing with another software (gvSIG or ArcGIS) it works fine.

Anyone knows what is happening?


texts: How to use Data Driven Pages to generate shapefiles instead of generating PDF?

texts: Is there any way to use Data Driven Pages to generate shapefiles instead of generating PDF? Because I need to export shapefiles generated with their layout in DWG format. If not, do you have an idea of the procedure that I can follow to achieve this task in a more effcient way?


texts: Script to “extract by mask” a list of raster files and list of mask

texts:
I'm new to Python so sorry for maybe trivial question esta, I have four Landsat images and each image with an area to extract by mask. How can I make this script? I've tried with this script doesn't work

import arcpy, os
from arcpy import env
path = 'F:\\Master\\Python\\1997\\data'
arcpy.env.overwriteOutput = True
arcpy.env.workspace = path
arcpy.CheckOutExtension("spatial")

Local variables:

1='1_23165_1997204.img'
2='2_23065_1997197.img'
3='3_23166_1997204.img'
4= '4_23066_1997197.img'
rec_1 = "rec_1"
rec_2 = "rec_2"
rec_3 = "rec_3"
rec_4 = "rec_4"
Zona1 = "F:\\Master\\Python\\output\\zona1"
Zona2 = "F:\\Master\\Python\\output\\zona2"
Zona3 = "F:\\Master\\Python\\output\\zona3"
Zona4 = "F:\\Master\\Python\\output\\zona4"

Process: Extract by Mask

arcpy.gp.ExtractByMask_sa(1, rec_1, Zona1)
arcpy.gp.ExtractByMask_sa(2, rec_2, Zona2)
arcpy.gp.ExtractByMask_sa(3, rec_3, Zona3)
arcpy.gp.ExtractByMask_sa(4, rec_4, Zona4)


texts:
I have a shapefile of population points, which I've categorized to show the difference in size for each point. However, I'm not satisfied with the range of values it has set out (in picture below).
Would there be a way to exclude values of 0 so that it starts at 0.1 up to 1 instead of 11? And so on for the rest of the points?


texts: Calling "fix geometries" in the QGIS python console?

texts:
I want to use the "fix geometries" tool in the Python console in QGIS.
However, I can't find anywhere, how to function is called and how to use it as code.
I tried
QgsProcessingOutputVectorLayer(input_path, output_path)

But I guess this is wrong, since it does not work. Any hints?


texts: Google Earth Engine: calculating and plotting monthly average NDVI for a region from MODIS data

texts:
I am trying to calculate monthly average NDVI for Ukraine the 2010-2020 period using MODIS data and then plot the values onto a graph (averaging all images in each individual month of each year separately across the full period). You can see the code I am using below. However, I keep getting the following error:

From checking the GEE documentation, I assumed this was an casting issue and have fixed those issues so I am not sure why the code is still failing.
//calculate average monthly NDVI for Ukraine
//get MODIS NDVI
var collection = ee.ImageCollection('MODIS/006/MOD13A1').select('NDVI');

//import country boundaries
var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');

//Subset Ukraine feature from countries.
var ukr = ee.Feature(
  countries
   .filter(ee.Filter.eq('country_na', 'Ukraine'))
    .first()
);

// Choose a reference period for your baseline
var reference = ee.ImageCollection(collection.filterDate('2010-01-01', '2020-12-31'))
  .sort('system:time_start');

//function to calculate monthly mean
//Get timestamp of first image
var time0 = ee.Image(reference.first().get('system:time_start'));

//Create an empty list (zeros) list to provide a structure into which the NDVI anomalies can be mapped
var first = ee.List([
  ee.Image(0).set('system:time_start', time0).rename('NDVI')
]);

//Use a function to iterate through months and calculate average NDVI
var months = ee.List.sequence(1, 12);
var byMonth = ee.ImageCollection.fromImages(
      months.map(function (m) {
        return ee.Image(reference.filter(ee.Filter.calendarRange(m, m, 'month')))
                    .select(1).mean()
                    .set('month', m);
}));

//Run the iteration
var average = ee.ImageCollection(ee.List(
  reference.iterate(byMonth, first)
));

//Define chart parameters
var chartParam = {
 title: 'Monthly average NDVI',
  hAxis: {title: 'Time'},
  vAxis: {title: 'Average NDVI'},
};

//Plot the chart
var chart = ui.Chart.image.seriesByRegion({
  imageCollection: average,
  regions: ukr,
  reducer: ee.Reducer.mean(),
  scale: 500,
  xProperty: 'system:time_start',
  seriesProperty: 'PROJECT'
});
//Print chart to console
print(chart.setOptions(chartParam));


texts:
I have City points and a Country polygons. Is it possible to make it so that labels of the cities are inside of their own country without converting to annotation?

I'm using ArcGIS 10.1 for Desktop.


texts: Cancel/pause a Geoprocessing Task without rolling back?

texts:
I have a Snap Geoprocessing task that is for about 57,000 points to lines. The issue is I need to restart the VM it is running on.

I have no idea how long it is going to run as it has been going for 48+ hours now. I do not want to cancel it and lose everything that has been done so far. Is there a way to cancel/pause it and save the edits that already occurred without it rolling back??

I'm fairly certain it is still running. I have a Fiber Manager extension and it is filling in a GUID as it is going, that is why it is taking long. It is also asking me to select from multiple objects once in a while, so it is not hanging, as of yet anyway.

It is an Enterprise Geodatabase 10.3.1, SQL Server 2008R2. I have tried to cancel it after a couple hours before and it rolled back.


texts: Extracting color values from polygons drawn in Google Earth Engine?

texts:
I have drawn polygons around Solar Panel units in California in google earth engine. I am interested in extracting the individual color values of the Red, Green and Blue bands within these polygons.

Anyone know how to do this?

I am thinking it is a relatively simple task but I am not sure where to look for these color values.


texts: Using a feature in a function in OpenLayers

texts:
I want to use in the get_pdf_data function the data loaded in the "imobile" vector so that I can perform operations with the "imobile" vector.
The code I'm trying to do this with is:
  var vectorSource = {
       teren: new ol.source.Vector(),
       constructii: new ol.source.Vector(),
       uat: new ol.source.Vector(),
       search: new ol.source.Vector()
 };

var imobile = new ol.layer.Vector({ source: vectorSource.teren, declutter: true });

      function get_pdf_data(idPdf, idPrimarie){
                console.log('DATA PDF');
                $('#container_respok').hide();
                $('.container_hide').show();
                $('.container_cauta').hide();

                $('#map').css('width', '50%');
                $('#map').css('float', 'right');

                var extent = ol.extent.createEmpty();

      imobile.getSource().forEachFeature(function(feature){
              feature.setId(feature.get('text'));

              var data = "";
              if($('#etc_tarla').is(":checked")){ data = "T:"+feature.get('tarla')+"  "; }
              if($('#etc_parcela').is(":checked")){ data = data+"P:"+feature.get('parcela')+"  "; }
              if($('#etc_suprafata').is(":checked")){ data = data+"S:"+feature.get('suprafata')+"MP.  "; }
              if($('#etc_proprietari').is(":checked")){ data = data+"P:"+feature.get('proprietari')+"  "; }
              if($('#etc_eterra').is(":checked")){ data = data+"C.F.:"+feature.get('e_terra')+"  "; }

         if(feature.get('text') == idPdf){
         console.log('GASITTTTTTTTTTTTTTTTT');

         var feature = imobile.getSource().getFeatureById(idPdf);
         //console.log(feature);

         selected_stile.getText().setText(data);
         feature.setStyle(selected_stile);
         ol.extent.extend(extent, feature.getGeometry().getExtent());
        }
      });

         map.getView().fit(extent, {size: map.getSize()});

     }

And the error received is the following:

Uncaught TypeError: imobile.getSource is not a function
at get_pdf_data (harta2.php: 3238)
at meniu_search (harta2.php: 3280)
at HTMLSelectElement.onchange (harta2.php: 1)

This is very useful to me because when I resize the map to 50% width automatically I have to look for a terrain in the respective map with idPdf and then to center the map on it


texts: Why can't I edit attribute table imported through 'text as layer'

texts:
I imported the coordinates from an Excel table. They show as single points as intended. Now I want to add columns to the attribute table, but I'm not able to edit it. I wanted to import the following columns from Excel as well, but that doesn't seem to work either.

Sorry for the noob-question, but I couldn't find anything relating the theme through google and such.

Ok, with pictures then:
I want to use that to add columns with information to these

You see, the small pensil thingy is greyed out, so I can't edit anything in this layer and I can't add any information to my points.

If I try Ale suggestion I get this error:


texts: Attempting to Run the Find Closest Facility Tool Within Python For Loop

texts:
While using ArcGIS Pro, I am trying to run the Find Closest Facility tool inside a Python For loop. I found some helper code in the tool information rollover in Pro that goes through these steps:

username = "<your user name>"
password = "<your password>"
cf_service = "https://www.arcgis.com".format(username, password)

#Add the geoprocessing service as a toolbox. Use an alias when importing
arcpy.ImportToolbox(cf_service, "agol")

#Set the variables to call the tool
facilities = r"C:/data/Inputs.gdb/Stores"
incidents = r"C:/data/Inputs.gdb/Customers"
output_routes = "C:/data/Results.gdb/Routes"
output_directions = "C:/data/output/Results.gdb/Directions"

#Call the tool
result = arcpy.agol.FindClosestFacilities(incidents, facilities, "Minutes", "", 2)

After entering my user/pass, I tried to import the Find Closest Facilities tool from AGOL and was presented with this error:

"OSError: The toolbox file http://logistics.arcgis.com/arcgis/services;World/ClosestFacilty;<username>;<password> was not found."

My user account is an Admin so I should have privileges to call this tool from AGOL to run inside of a Python script. Perhaps, I'm calling it incorrectly? I tried using the URL above, as well as, our AGOL portal URL.

What I'm ultimately trying to do is to just run the "FindClosestFacilties" tool in a Python script, so if anyone has an alternative to the above Esri code, that would be helpful.


texts: QGIS: interpolate difference in Z-value on enbankments

texts:
I have two layers:

1 polygon layer with a theoretical profile of a harbour.
1 pointlayer with output from my multibeam echousounder.

I made a new pointlayer with a Z-difference between the two layers.
The output is a value above or under the theoretical profile by comparing the Z from a point to the Z from the theoretical profile using [geomwithin]

this is the result chart:

The red spots are above the theoretical profile, the green spots under the theoretical profile.

Here is my problem.

The blue lines are the outlines of the polygons. I'm good with the polygons that have only 1 value (the flat areas in the theoretical profile), but I have a problem with sloped embankments.

The provider of the theoretical profile put a min and max z value in the data. But it is not specified in a way that I can use it.

My goal:

Make a chart that also calculates the embankments (the slopes) Z-difference, instead of the current method where the slopes are only compared to 1 value.

I would like to calculate the green areas.

Am I in the right direction in thinking:

I need to add the right Z-value to the vertices of the polygon of the theoretical profile
I need a method to interpolate between the vertices with a minimal and maximum Z-value of the theoretical profile.
I need a method to compare the measured Z-value and the interpolated Z of the slopes.

I would like some advice on the validity of my method. and some guidance in how to interpolate or make a 3D version of the theoretical profile.


texts: Changing the CRS of all shapefiles in QGIS

texts:
I've just realised I've been using the wrong CRS for an analysis.

As I've done a fair bit of analysis (distance to nearest hub, transformation of lines into points etc), I was wondering if by setting the new CRS I need, it would update all the distances inside the shapefiles etc or do I have to run everything again?

For example, I have the distance of each city from a border, set in ESPG::26985, and I'd like to change it in ESRI::53027. I have this for 25 countries...


texts: WFS-T returns 'Success' but doesn't delete feature

texts:
I have a layer on my local GeoServer. I have published it, now I can successfully access it with WFS operations (DescribeFeatureType, GetFeature etc.). I have also enabled the transactions on WFS and give the necessary write permission to all the roles (including the anonymous roles). I am trying to delete a feature from my layer by using WFS-T.

Here is my URL (Aydin1 is the workspace name):
http://localhost:8080/geoserver/Aydin1/ows

Here is my WFS-T XML which I use in my POST request:

<?xml version="1.0" ?>
   <wfs:Transaction
   service="WFS"
   version="1.0.0"
   xmlns:wfs="http://www.opengis.net/wfs">
   <wfs:Delete typeName="Aydin1:PointLayer">
        <Filter>
            <PropertyIsEqualTo>
                <ValueReference>details</ValueReference>
                <Literal>DetailsPoint1</Literal>
            </PropertyIsEqualTo>
        </Filter>
   </wfs:Delete>
</wfs:Transaction>

And here is the answer I get from the server:

<?xml version="1.0" encoding="UTF-8"?>
<wfs:WFS_TransactionResponse version="1.0.0"
    xmlns:wfs="http://www.opengis.net/wfs"
    xmlns:ogc="http://www.opengis.net/ogc"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    xsi:schemaLocation="http://www.opengis.net/wfs http://localhost:8080/geoserver/schemas/wfs/1.0.0/WFS-transaction.xsd">
    <wfs:InsertResult>
        <ogc:FeatureId fid="none"/>
    </wfs:InsertResult>
    <wfs:TransactionResult>
        <wfs:Status>
            <wfs:SUCCESS/>
        </wfs:Status>
    </wfs:TransactionResult>
</wfs:WFS_TransactionResponse>

So, as you see, my transaction is successful, I don't get any error message. However, when I check my layer on GeoServer, I still see my feature, it hasn't been deleted. I know that my Filter structure is just fine; the exact filter is working well when I use GetFeature WFS operation.

I have published two different layers. One of them is based on shapefiles, the other one comes from PostGIS. They act the same way, that is, I get the SUCCESS response message but nothing happens.

What could be the problem?


texts: I'm trying to find an equivalent in FME to the Majority Filter in ArcGIS (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/majority-filter.htm), I've had no success this far... Help, anyone?


texts: Not able to save edits to vector layer in QGIS

texts: I have a vector layer with about 200 features. So far I have been able to add polygons to this layer just fine, now it is giving me trouble. When I draw the polygon and try to save edits, the polygon disappears. In a few cases the polygon is saved successfully though it is rare. Then occassionaly the polygon disappears but the attributes that are assigned do appear in the attribute table. Any suggestions how to get through this?


texts:
I am recently working on a project where i want to show a skiarea in a perspective view, like you know it from maps of skiresorts where the ski lifts are displayed. Are there any libraries where i can achieve this or must i do this by myself?

I want to create something like this, but autogenerated (could be a picture). I want to display it on Android or iOS

I researched a lot and I think there is not such a thing, but you can prove me wrong, hopefully.


texts: When the table is too big for it to get exported to excel in arcgis ModelBuilder

texts:
My model output got way too many rows to be exported to excel in arcgis ModelBuilder. It would work if I could use the new xlsx file format but it seems the toolbox can only do xls. After a little research I found that cvs might solve my problem, however when I attempted to use it I got the following error message;

Start Time: Tue May 24 08:37:19 2016
Running script ExportXYv...
Writing Results to ASCII File....
Failed script ExportXYv...

Traceback (most recent call last):
  File "c:\program files (x86)\arcgis\desktop10.3\ArcToolbox\Scripts\ExportXYV.py", line 159, in <module>
    export = setupXYV()
  File "c:\program files (x86)\arcgis\desktop10.3\ArcToolbox\Scripts\ExportXYV.py", line 50, in setupXYV
    outFieldNames = outFieldNames)
  File "c:\program files (x86)\arcgis\desktop10.3\ArcToolbox\Scripts\ExportXYV.py", line 139, in exportXYV
    UTILS.writeText(fo, "{0}\n".format(outRow))
UnicodeEncodeError: 'ascii' codec can't encode character u'\xc5' in position 624: ordinal not in range(128)

Failed to execute (Export Feature Attribute to ASCII).
Failed at Tue May 24 08:37:20 2016 (Elapsed Time: 0,79 seconds)

Does anyone know how I can solve this problem?

When I attempt to use the "table to excel" tool I get the following error:

"Executing (Table To Excel): TableToExcel C:\Users\fmromss\Documents\ArcGIS\Default.gdb\ro_arealressurs_f_fkb_ny_326_union C:\Users\fmromss\Documents\ArcGIS\scratch\ro_arealressurs_f_fkb_ny_326_union_TableToExcel.xls NAME CODE
Start Time: Tue May 24 08:59:14 2016
Running script TableToExcel...
ERROR 001531: Input table exceeds the 65535 rows limit of the .xls file format.
Completed script TableToExcel...
Failed to execute (Table To Excel).
Failed at Tue May 24 08:59:14 2016 (Elapsed Time: 0,82 seconds)"

As it takes a lot of time to scroll in the table Im not completely sure about how many rows we are talking about here but I know its over 180 000


texts: QGIS buffer of lines results in "random" vertex order in buffer polygon

texts:
I tried this topic already on QGIS dev mailing list:
I have a fairly complex geometry-(and mind-)bending workflow where I'd ideally depend on the fact that buffering a line (regular buffer, flat cap, miter join with almost infinite limit) would result in a reproducible pattern for the vertex order in the resulting buffer polygon (s. [1] for details).
I found that in 99% of the cases the first vertex of the buffered polygon will be on the left side of the original line feature's second vertex (which is great). That's confusing, not sure how else to put it, but see this image for viz, showing the vertex order of the original line features and the resulting polygon features:

However, on a few occasions this assumption doesn't seem to hold and the first vertex of the buffer polygon will be anywhere with respect to the original line feature's first vertex and does not follow any pattern at all, s. [2] for a WKT example in EPSG:25832 and buffered with 0.3 m (it's too long for an image). On this feature, vertex 0 of the polygon feature is on the left side of vertex 14 of the line feature, not next to vertex 2 as for all other line features I buffered.
My question is: is that expected behavior that buffer feature vertices can start at a sorta random place with respect to the line feature's first vertex? Or can this be considered a minor bug?
[1] I want only one side of the buffered geometry and need to somehow rotate the buffer vertex list so that the first vertex is related to the first vertex of the line feature. I can't use singleSidedBuffer() because that is much more likely to produce weird artifacts on sharp kinks in the line geometry
[2] The WKT of a feature where it's not working as I'd expect (in EPSG:25832, buffered with 30 centimeters, flat cap, miter join and 100m miter limit): LineString (567128.05779425 5674382.50117657, 567124.28299077 5674385.77951107, 567121.44049412 5674387.67067814, 567113.39428435 5674393.02397452, 567089.71334818 5674408.77935169, 567075.32187873 5674418.35427, 567048.09794292 5674433.92388853, 567047.26883028 5674434.39806565, 567046.22311791 5674435.0838905, 567043.6035418 5674436.80192543, 567006.05562972 5674453.08412973, 566973.94524544 5674467.01295346, 566941.83486115 5674480.9417772, 566930.00051817 5674486.07016386, 566929.97872742 5674486.07960683, 566929.33462488 5674485.86262721, 566904.04808487 5674477.3443191, 566871.1140052 5674462.85053673, 566838.17992554 5674448.35675436, 566797.46760646 5674431.84405364, 566756.75528738 5674415.33135292, 566716.04296829 5674398.8186522, 566675.33064921 5674382.30595149)


texts: How to read/write from/to QGIS project file using Python?

texts: I'd like to save the project-dependent configuration of my Python plugin inside the QGIS project file. Is this possible and what are best practices to do so?


texts: Transparent background after warp operation that rotate the image

texts:
I am trying to reproject a georeference image using gdalwarp (and after directly in c++). The reprojection implies a rotation, so the resulting image has a black background where I want to have a transparent background.

The best I could find was to run something like that:

gdalwarp -srcnodata 255 -t_srs EPSG:4326 test.kap test.tif

But this also remove black pixels in the middle of the map.

I have also tried many other things and options without any luck.

Is there a way to achieve this? If not is there a way to calculate the raster "footprint" so I can use it as a mask to hide pixels with no data?

I have tried many options including:

gdalwarp -dstalpha -t_srs EPSG:4326 test.kap test.tif

The result is that the background is transparent but the map itself is just plain green (as I can see in QGIS).


texts:
Im sure there is a simple answer to this but its eluding me. When saving outputs/QGIS projects the file navigator interface opens up and I can go through my files until I find the place I want to save the work/project. However, my machine has two disks, a 256GB SSD and a 2TB HD. The QGIS installation is on the SSD but I save all of my work on my HD because of its larger capacity. However, I can't change volumes when I want to save my work. It will allow me to navigate the SSD but I can't see where I can change the volume Im saving the work to. I could save it in a temp folder and move it later, but its a bit of a pain. When I click on the volumes under 'favourites' or 'recently used', only the SSD is listed.

Is there a way to change the directory to the volume I want?

I have the latest QGIS (pisa) on a Mac OS


texts: How to blend borders of multiple raster files so that it doesn't leave any discrepancies while mosaicing them together?

texts: I am using landsat data to create an FCC over a large basin area. I combined the raster files to get the entire area of the basin as one raster. But as the images were taken on different days, they don't merge together as one but look as if multiple images have been put together. Is there any way to blend the individual images into a whole ?


texts:
I've encountered an issue: when I try to execute in Spatialite_GUI query that contains WITH clause I get SQL error: near "with": syntax error. Also WITH unlike SELECT and other keywords is not highlighted inside SQL-querry window. The same query in SQLite Manager works just fine. Found no treads on this issue. Is this my local problem?

Spatialite_GUI 1.7.1-1.5, SQLite Manager 0.8.1, OpenSUSE 13.1.

EDIT:

Actually the issue is deeper. DB Manager in QGIS 2.4 also throws the same error for the WITH clause (WHITH in sql window is highlighted in this case).

  General info

  SQLite version:
  3.7.17

  SpatiaLite

  Library:
  4.1.1  GEOS:
  3.4.2-CAPI-1.8.2 r3921  Proj:  Rel. 4.8.0, 6 March 2012


texts: QGIS Server: SLD and SLD_BODY parameters in GetMap request

texts:
Does Qgis server support SLD and SLD_BODY parameters in GetMap request?

I'm trying to give a valid .sld to qgis server but it returns empty image in response.