texts: Solve Esri Toolbox produces empty Shapefiles (make route layer)

texts:
I am making a Toolbox for finding the best Pizzaroute between 7 Pizzastops.

We used the make.routelayer from ArcPy. Our Pizzaroutetime expressed not our expression, rather all Firestations (station 1-49). Our Outputeroute shows a route between all Firestations.

How do we get the right output?

import arcpy
from arcpy import env
import time
import os
arcpy.env.workspace = "d:\Users\XXX\Desktop\NetworkAnalyst\Tutorial\SanFrancisco.gdb"
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("network")

NetzwerkWork = arcpy.GetParameterAsText(0)
  if NetzwerkWork == '':
     NetzwerkWork = 'd:\Users\XXX\Desktop\NetworkAnalyst\Tutorial\SanFrancisco.gdb\FireStations.shp'
arcpy.env.workspace = NetzwerkWork

networkDataset = "d:\Users\XXX\Desktop\NetworkAnalyst\Tutorial\SanFrancisco.gdb\Transportation\Streets_ND"
stops = "d:\Users\XXX\Desktop\NetworkAnalyst\Tutorial\FireStations.shp"
arcpy.env.overwriteOutput = True
stops_select = "d:\Users\XXX\Desktop\NetworkAnalyst\Tutorial\SanFrancisco.gdb\Pizza_Time"
stops_select2 = "d:\Users\XXX\Desktop\NetworkAnalyst\Tutorial\SanFrancisco.gdb\Pizza_Distance"

Adress_locator = arcpy.GetParameterAsText(1)

Punkt_1 = arcpy.GetParameterAsText(2)

Punkt_2 = arcpy.GetParameterAsText(3)

Punkt_3 = arcpy.GetParameterAsText(4)

Punkt_4 = arcpy.GetParameterAsText(5)

Punkt_5 = arcpy.GetParameterAsText(6)

Punkt_6 = arcpy.GetParameterAsText(7)

Expression ="NAME = 'Station 1'" + " ""OR" +" " "NAME = 'Station 2'"+ " ""OR" +" " "NAME = 'Station 3'"+ " ""OR" +" " "NAME = 'Punkt_3'"+ " ""OR" +" " "NAME = 'Punkt_4'"+ " ""OR" +" " "NAME = 'Punkt_5'"+ " ""OR" +" " "NAME = 'Punkt_6'"
Pizza = arcpy.Select_analysis(stops, "Pizza_Time", Expression)
Pizzaroute_Time = arcpy.MakeRouteLayer_na(networkDataset, Pizza, "TravelTime", "USE_INPUT_ORDER", "PRESERVE_BOTH", "NO_TIMEWINDOWS", "", "ALLOW_UTURNS", "Oneway;RestrictedTurns", "USE_HIERARCHY", "", "TRUE_LINES_WITH_MEASURES", "").getOutput(0)

Timeclass = arcpy.na.GetNAClassNames(Pizzaroute_Time)

TimeSublayer = arcpy.mapping.ListLayers(Pizzaroute_Time, Timeclass["Routes"])[0]

arcpy.na.AddLocations(Pizzaroute_Time, Timeclass["Stops"], stops)

arcpy.na.Solve(Pizzaroute_Time)

arcpy.management.CopyFeatures(TimeSublayer, "Pizza")

arcpy.AddMessage("Completed")

Additional information after receiving answer...

The SQL expression works fine.
Our Actually Expression looks like this:

Expression = "NAME = 'Station 1'" + " ""OR" +" " "NAME = 'Station 2'"+ " ""OR" +" " "NAME = 'Station 3'"

And we get an Output including only these 3 Points (Pizza_Time).
Our Problem is that we can not do a route between these 3 Points from the Expression, Because the MakeRouteLayer Output ("Pizza") gives us a route between all the Points of the Input-Shapefile (FireStations).


texts:
I'm testing the Script Runner plugin to run scripts in QGIS (I'm new with programming with Python in QGIS).
For some reason I can't get a script with arguments to work.
It throws this error with any script I execute (the sample ones too):

     Traceback (most recent call last):
       File "C:/Users/ag/.qgis2/python/plugins\scriptrunner\scriptrunner.py", line 548, in run_script
         print "----------%s----------" % datetime.datetime.now()
       File "C:/Users/ag/.qgis2/python/plugins\scriptrunner\stdout_textwidget.py", line 55, in write
         sys.__stdout__.flush()
     IOError: [Errno 9] Bad file descriptor

I've been checking both scriptrunner.py and stdout_textwidget.py, but I can't see what's the problem.


texts: Trying to identify the first upstream point from each point in a feature class (geometric network)

texts:
I'm new to geometric networks, but someone suggested it on a different question and it seems like it might have the capability I'm looking for. I see that within a geometric network you can trace upstream, and find all points that are upstream from the points you select. I need to find the one point from one feature class (electrical power sources) that is the first upstream from each point in a separate feature class (houses). So each house point will have an associated power source point, so that I may add an attribute to the house points populated with the ID number of the electrical source that feeds it. I'm going for an outcome similar to a spatial join, but instead of intersecting or being within a distance, the join features are the first points upstream from the target points. Right now I only see how to place flags manually (I need several thousand) and the results of the trace are from the same feature class, and select everything upstream, not just the first.
Is it possible to trace upstream for every point in a feature class, and have the results be from a separate point feature class? I am using ArcMap 10.6 basic license


texts: I want to extract the Z values of polygons from a GML file with a python script. I can't use the field calculator, since you can't edit GML files in QGIS. Ideally, I want to extract the Z values from each node and then take the mean of these values and store it in a list in Python.


texts: Repeatable analysis for new maps on ArcGIS Online

texts:
I need to find a way to have a repeatable analysis in ArcGIS Online that can be applied to multiple layers uploaded separately.
Each time there is a new run on our software, the client will have the option to send the results to ArcGIS Online. These will appear as two layers within the same file and will always have the same attribute columns. The user will always want to see the same visualization of these results but it takes some time to set up- how can we just re-apply the same analysis between files?
Very new to ArcGIS.


texts: How to buffer a pair of longitude and latitude coordinates in R?

texts:
I have a pair of longitude and latitude coordinates, and I want to plot a 500 meter buffer around it using leaflet or mapview. st_buffer, however, does not correctly buffer longitude/latitude data.

Initial code:

library(rgeos)
library(rgdal)
library(sf)

df <- data.frame(lon = -121.9552, lat = 37.35411)
coordinates(df) <- c("lon", "lat")
proj4string(df) <- CRS("+init=epsg:4326")
df_sf <- st_as_sf(df)
#....
# Something goes here, but I don't know what?
#...
st_buffer(df_sf, 500)

Warning message:
In st_buffer.sfc(st_geometry(x), dist, nQuadSegs) :
st_buffer does not correctly buffer longitude/latitude data

What's the missing piece?


texts: France, georeferencing with postal codes not working in CartoDB

texts: CartoDB is not georeferencing many of my locations in France. I am using postal codes, and specifying all of them are in France. CartoDB is able to georeference only about half of them. Any ideas?


texts: How do I convert the LST values on the MODIS LST Image to degree celsius

texts: I have downloaded MODIS 8 day average LST images, but the values are way too high. How do I convert the pixel values to degree celsius? Someone said they multiplied the data with 0,02 and then subtracted 273 to get them as Celsius. is that the correct method?


texts: I need to recreate maps like on picture in QGIS 3. They were created in ArcGIS 9.3. Have no idea how to do this, found no guides. The main problem are roads and rivers, because they are not specific dots, but line. Maybe some plugins? I am pretty new to GIS systems.


texts: How do I convert WGS84 lat,long points from degrees to meters in Python?

texts:
I am facing a problem in the creation of heat map based on the lat-long points. All the points are closely packed into a group like squares, so I would like to convert the lat,long points to meters so that I can plot it over the graph to see the variation in the heatmap.

Current my points are in WGS84 (4326) format.
How can I convert projection coordinates into the meters so that I can plot the heat map?


texts: Make only one attribute editable in a webmap

texts:
Can someone tell me (spoon feed me) how to make only one attribute editable in ESRI online. Everything I have read and followed makes all the attributes editable.

I really need the public user to enter a date in the contractor start date attribute field, it would be best if they were not allowed to delete anything.


texts: ogr2ogr of NTF_Lambert_II_etendu is off by about 30 meters

texts:
I am converting RPG files in NTF_Lambert_II_etendu to wgs_84 using:

ogr2ogr -t_srs EPSG:4326 ilot_2008_027_wgs84_2.shp ilot_2008_027.shp

but the result is off by about 30 meters. Surprisingly, the same conversion using Qgis works fine. Could anybody tell me how to correct that or how to investigate the exact transformation qgis is doing?

The proj file of the intial file is :

PROJCS["NTF_Lambert_II_etendu",GEOGCS["GCS_NTF",DATUM["D_NTF",SPHEROID["Clarke_1880_IGN",6378249.2,293.46602]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",600000.0],PARAMETER["False_Northing",2200000.0],PARAMETER["Central_Meridian",2.3372291667],PARAMETER["Standard_Parallel_1",45.8989188889],PARAMETER["Standard_Parallel_2",47.6960144444],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",46.8],UNIT["Meter",1.0]]


texts: I am working with Sentinel 2 bands in Qgis 2.18. When I load them to Qgis' TOC, they appear in black and in the layers panel they do not display the minimum and maximum number of pixels. I would like to solve this problem.


texts: Why C# Server Object Extensions are faster

texts:
From my experience C# SOE are the fastest. It is also discussed here.

ArcGIS Server is developed in Java, how come .NET is the fastest?


texts:
I have PostGIS 2.0.1 and PgRouting latest version. I'm using osm2pgrouting to import OpenStreetMap data everything goes OK the trouble is that I always get source and target columns without data and when I run SELECT assign_vertex_id('ways', 0.00001, 'the_geom', 'gid'); by hand on PGSQL I get the following error:

  ERROR:  syntax error at or near "vertices_tmp" LINE 1: SELECT
  ST_addGeometryColumn(''vertices_tmp'', ''the_geom'', ...
                                        ^

  ********** Error **********

  ERROR: syntax error at or near "vertices_tmp" SQL state: 42601
  Character: 31

After hours of reading on the net I've realized that the function assign_vertex_id() had to change somehow because it references other functions like StartPoint that it's called ST_StartPoint and I've changed that but I still get errors.

Where can I get the .sql file with the proper SQL to make the upgrade and postgis/pgrouting can work together?


texts: How to omit NA values in the final plot when plotting a SpatialPolygonsDataFrame in R?

texts:
I'm creating a map of Pakistan based on availability of satisfactory water supply at a district level.

Since it is difficult to collect reliable data for all districts, I would like for the final map to show whitespace where the column Condition is unavailable. My SPDF is named pak.

tm_shape(pak) + tm_polygons("Condition", palette = "-Reds")

Above is the command I am using presently to get the following output. I would essentially like for the grey space to not be plotted at all (not even the borders).


texts: Retaining input information with raster overlay

texts:
I'm working in ArcGIS (version 10.7) on an unweighted overlay analysis of 15 green space functions that I would like to tally but also still retain the information for which functions are represented within different landscape units, as well as the total number of functions represented within each landscape unit.
I've rasterized and classified each of the 15 functions into a discrete scoring system (0= functions do not exist, 1= function exists) and used raster calculator to produce a raster with the total value/function counts (0-15), but that of course just gives me the total count and not the information about which functions are represented.
Does anyone know how I can most efficiently join the landscape units and information about functions represented?


texts:
Using QGIS 2.12.1-Lyon on OSX Yosemite, I am following this tutorial.

I cannot add any polygon on the map : I can draw the polygon, but as soon as I right-click, the polygon disepear, even if the Features Attributes windows is prompted.

Editing the feature (Polygon is stretched in order to reach the screenshot tool)

After right-clicking

Am I doing something wrong or is there a bug ?


texts: ArcGIS Online Feature Style Using Multiple Attributes

texts: I am configuring style for a line feature in ArcGIS Online. I used one string field to define the line style and want to use another string field to define the color. However the second attribute is not in the drop-down list. Is there any requirement for attribute field to be used as a second option for styling? How can I do this?


texts: Draw label outside polygon, if it not fit's inside (MultiPolygon, QGIS)

texts:
I would like to label polygons (parts of multipolygons) inside IF the label fits completely inside, else the label should be outside with a call out.
I tryed that, as you can see in the screenshot:
I was using a expression, to force the drawing of the label 4mm distance to anchor point if the area is smaller then 25ha. Works with #128 but not with #106 - so this solution is not working. (Even if that would work, it is not a good solution since it's scale depended!)

The task seems simple, but I found no solution: Draw label outside, if it not fit's inside. If it's outside, draw it with call out to the objects interior. (#128 is the expected result)


texts: ArcGIS Pro performance issue when using Text in Layout?

texts:
A few weeks ago I began to notice that certain map projects I was working on in ArcGIS Pro were running extremely slow when I tried making adjustments to the map Layout, specifically moving around text elements. The confusing thing for me is that the functions you would expect to make the system lag, like geoprocessing or editing, are running fine; it doesn't appear the content of the maps has anything to do with the problem. The issue only appears when I try to add or move around text elements in the Layout pane.
Has anyone experienced anything similar and have a fix for this?


texts: How can I find the nearest points (plural) on a line segment using PostGIS

texts:
I'm locating the nearest point on a line segment with PostGIS and ST_ClosestPoint, I can select the point on each of these line segments which is closest to the starting point. Which looks something like this:

But! I cannot figure out how to get the second closest point, or the third closest point on a given segment.

What I would like is more like this:

I'd like a range of nearby points on the line segments.

Any thoughts about how this might be accomplished?


texts: Merge point features within a given radius

texts: I have a file in QGIS with 5000 points. Some of the points are really close to each other (5m). How can I cluster the point features to one point feature if they are within 5m? The new points should be centroids from the given points.


texts:
I have created a project where I need to add and remove multiple layers in OpenLayers. The layers are from shapefiles that are available on a node server. I loop trough a list of all the names of the files and add them one by one to the map.
  setShapeFile(shapefiles) {
    this.clearAllBoundriesLayers();
    this.resetShapefileColors();
    shapefiles.forEach((shapefile) => {
      let vectorSource = new VectorSource({
        format: new GeoJSON(),
        loader: (extent, resolution, projection) => {
          var url = "http://localhost:3030/api/getShapefile/" + shapefile.file;
          var xhr = new XMLHttpRequest();
          let color = this.getshapefileColor();
          xhr.open("GET", url);
          var onError = function () {
            vectorSource.removeLoadedExtent(extent);
          };
          xhr.onerror = onError;
          xhr.onload = () => {
            if (xhr.status === 200) {
              var features = vectorSource
                .getFormat()
                .readFeatures(xhr.responseText);
              features.forEach((feature) => {
                feature.setId(feature.get("OBJ_ID"));
                feature.setStyle(
                  new Style({
                    fill: new Fill({
                      color: color,
                    }),
                    stroke: new Stroke({
                      color: "#000000",
                      width: 0.5,
                    }),
                  })
                );
              });
              vectorSource.addFeatures(features);
            } else {
              onError();
            }
          };
          xhr.send();
        },
        strategy: bboxStrategy,
      });
      let vector = new Vector({
        //minZoom: 13,
        source: vectorSource,
      });
      vector.set("name", "plotShapefileLayer" + shapefile.type);
      this.map.addLayer(vector);
    });
    console.log("-----------------------");
    this.map.getLayers().forEach((layer) => {
      console.log(layer);
      console.log(layer.get("name"));
    });
  }

In the console log I do at the end all of my layers are visible with the correct names.
But when I try to remove the layers only the first layer gets removed and the other ones stay visible. The layer objects become unidentified. So I was thinking that by removing the first layer all the other ones there reference to the layer gets removed but the layer self stays on the map. I have tested if the method to remove the layers doesn't get accidentally called, but this doesn't seem to be the case. Here is the code where I remove the layers:
clearAllBoundriesLayers() {
    console.log("-----------------------");
    this.map.getLayers().forEach((layer) => {
      console.log(layer);
      if (layer !== undefined) {
        console.log(layer.get("name"));
        if (layer.get("name") !== undefined) {
          if (
            layer.get("name") === "plotBoundriesLayer" ||
            layer.get("name") === "plotUserBoundriesLayer" ||
            layer.get("name") === "plotUserBoundriesLayerIcons" ||
            layer.get("name").includes("plotShapefileLayer")
          ) {
            layer.getSource().clear();
            this.map.removeLayer(layer);
          }
        }
      } else {
        this.map.removeLayer(layer);
      }
    });
    if (this.select) {
      this.select.getFeatures().clear();
    }
  }


texts: Using pyproj inverse transform with GeoPandas and GroupBy + SortBy constraints

texts:
I am using pyproj inverse transform to add azimuth and distance "info" to a grouped and ordered geodataframe. I asked a similar question earlier here, but I didn't realize grouping was a constraint, and I cannot figure out how to incorporate GroupBy into that answer. (The code in that answer is very compact and a bit hard to understand for a novice).
In the example below, I want to perform pyproj inverse transform in 'myid' groupings. In other words, I don't care about the distance and azimuth between point "1-3" and point "2-1"; that is to say, the distances and azimuths for every "X-1" point should be nan. How can I achieve this?

%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
import pyproj
from pyproj import CRS

#### Build example GeoDataFrame for easily reproducible example ####
myUTMepsg = 32611
myid = [1, 1, 1, 2, 2]
myorder = [3, 2, 1, 2, 1]
lat = [36.42, 36.4, 36.32, 36.28, 36.08]
long = [-118.11, -118.12, -118.07, -117.95, -117.95]
df = pd.DataFrame(list(zip(myid, myorder, lat, long)),
                  columns =['myid', 'myorder', 'lat', 'long'])
df.sort_values(by=['myid', 'myorder'], inplace=True)
df.reset_index(drop=True, inplace=True)
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['long'], df['lat']))
gdf_pt = gdf_pt.set_crs(epsg=4326)
display(gdf_pt.style.hide_index())
ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);
#### End

#### this is not necessary, but I like to add labels to the
#### points on the plot so that I can have a visual check
for x,y,z1,z2 in zip(gdf_pt.long, gdf_pt.lat, gdf_pt.myid, gdf_pt.myorder):
    label = str(int(z1)) + '-' + str(int(z2))
    plt.annotate(label, (x,y), textcoords = 'offset points',
                 xytext = (12,-5), ha = 'center')
#### End

#### Analysis & Problem Area
geod = CRS.from_epsg(myUTMepsg).get_geod()
for i, r in gdf_pt.iloc[1:].iterrows():
    myinfo = geod.inv(gdf_pt.long[i], gdf_pt.lat[i], gdf_pt.long[i-1], gdf_pt.lat[i-1])
    gdf_pt.loc[i, 'az_fwd'] = myinfo[0]
    gdf_pt.loc[i, 'az_back'] = myinfo[1]
    gdf_pt.loc[i, 'dist_meters'] = myinfo[2]
    gdf_pt.loc[i, 'bearing_degrees'] = max(myinfo[1], myinfo[0])

display(gdf_pt.style.hide_index())

I tried breaking out a working geodataframe inside the loop, but I ran into trouble trying to combine the working and final geodataframe; I kept getting duplicate entries (yellow below; red arrow is the one I want to keep), even when specifying keys, using Pandas merge, concatenate, and join.
#### Analysis & Problem Area
geod = CRS.from_epsg(myUTMepsg).get_geod()
for idval in gdf_pt.myid.unique():
    tempgdf = gdf_pt[gdf_pt.myid == idval].copy()
    tempgdf.sort_values(by=['myid', 'myorder'], inplace=True)
    tempgdf.reset_index(drop=True, inplace=True)
    for i, r in tempgdf.iloc[1:].iterrows():
        myinfo = geod.inv(tempgdf.long[i], tempgdf.lat[i], tempgdf.long[i-1], tempgdf.lat[i-1])
        tempgdf.loc[i, 'az_fwd'] = myinfo[0]
        tempgdf.loc[i, 'az_back'] = myinfo[1]
        tempgdf.loc[i, 'dist_meters'] = myinfo[2]
        tempgdf.loc[i, 'bearing_degrees'] = max(myinfo[1], myinfo[0])
    gdf_pt = pd.merge(gdf_pt, tempgdf, how='outer') ## I tried various key controls here without success
    del tempgdf

display(gdf_pt.style.hide_index())

Using Windows 10, conda 4.8.2, Python 3.8.3


texts:
I want to create a device that should transmit its location to a server using anyother medium other than GSM. Bascially I need a tracking device which i could give to my worker and he goes to a forest or mine he should be traceable without using any other medium like GSM. I want the satellite to monitor him.

Please tell me way


texts:
I am fairly new to GIS.

I have a line of a survey track that is approx. 380km long. I can correctly measure this using $length, and with the measure tool. I am using a custom transverse mercator custom CRS of my survey area.

When I try to place regular points along the line, using Points along geometry or QChainage, the algorithms seem to think my line is only a few metres long. If I space points every 'one metre', I will get 5 points, but if I try to space points every km, I will only get one point at the start of the line.

How can I get points every km?


texts: Is a plugin for QGIS3 with the user interface defined by code acceptable for the repository?

texts:
I am working on a QGis3x plugin, I have created a minimal Plugin, the user interface (dialog) I have defined by code. Without using Qt Designer.
If I would like to upload this plugin to the QGIS repository, For free use of the public, would it be acceptable?
If not, how can I implement QTranslator in this case?


texts: Calculating area between two SpatialLines

texts:
I'm trying to calculate the area between two lines. The two lines have the same start and end location, but the number of vertices for each line can vary.
I've tried to create a polygon from the coordinates of the vertices from the two line and then calculate the area by using rgeos::gArea but the resultant SpatialPolygons (sps) is incorrect.
x1 <- c(1,5,4,50)
y1 <- c(1,3,4,50)
line1 <- sp::SpatialLines(list(sp::Lines(sp::Line(cbind(x1,y1)), ID='a')))
x2 <- c(1,5,5,50)
y2 <- c(1,4,6,50)
line2 <- sp::SpatialLines(list(sp::Lines(sp::Line(cbind(x2,y2)), ID='b')))

plot(line1)
plot(line2, add = T, col = "red")

line_pts <- methods::as(line1, "SpatialPoints")
comparison_pts <- methods::as(line2, "SpatialPoints")

line_coords <- sp::coordinates(line_pts)
comparison_coords <- sp::coordinates(comparison_pts)

start <- line_coords[1, ]
end <- line_coords[base::nrow(line_coords), ]

coords <- base::rbind(line_coords, comparison_coords)

p = sp::Polygon(coords)
ps = sp::Polygons(list(p), 1)
sps = sp::SpatialPolygons(list(ps))

plot(sps, add = T, col  ="black")

Any ideas would be appreciated.


texts: Export MultiPolygon from MultiSurface layer in QGIS

texts:
I have a layer in QGIS that is a Multisurface layer currently stored on a PostGIS server I want to add this layer as a MultiPolygon layer to PostGIS using QGIS.
Can I specify this somehow when I save/export to QGIS? The only workaround I found is exporting the layer to a shape file loading it and exporting it then to PostGIS.


texts: QGIS Extract by Expression with multiple polylines in Processing Modeler

texts:
I'm trying to get the matching features from an intersecting vector layer. I can manually open Extract by expression tool and put the expression like intersects($geometry, geometry(get_feature('Liness', 'id', '1'))) where 'Liness' is one of the polylines and I'm giving the other polyline as input. I would like to implement it using QGIS Processing Modeler.

Since in extract by expression can only give one input, I'm unable to get the desired output. Can anybody suggest a better method? I want to use the matching features to other further steps.

Or can I implement it using pyqgis? Is there any code sample to refer?

EDIT 1: As suggested by Gabriel De Luca, tried "Extract by Location". But it's selecting my input polyline in random order.

Also, I would like to have each polyline with "id" of "Liness"(whose id changes based on the lines I plot for selection) to have a separate section. Why because, I would like to use its output to autoincrement function, where I can increase each of the input lines (whose selection came from different "id") values from zeros to a finite value based on the "Liness" intersection in order.


texts: Downloading MODIS data using Python with help of shapefiles

texts:
I want to download MODIS data at 250m (2 bands) and 500m (5 bands) for arond 2700 regions. I have the shapefiles for each of the region.

How can I download the raw MODIS data (surface reflectance) using Python (no GUI GIS software to be involved) using shapefiles (Esri format)?

I have also calculated the convex hull for each shapefile region, in case that is needed.

The code below uses latitude and longitude to cover the region - but this is not neat. It can miss tiles. besides this doesn't use the shapefiles.

import ee
import time
import sys
import numpy as np
import pandas as pd
import itertools
import os
import urllib

ee.Initialize()

def export_oneimage(img,folder,name,region,scale,crs):
  task = ee.batch.Export.image(img, name, {
      'driveFolder':folder,
      'driveFileNamePrefix':name,
      'region': region,
      'scale':scale,
      'crs':crs
  })
  task.start()
  while task.status()['state'] == 'RUNNING':
    print 'Running...'
    # Perhaps task.cancel() at some point.
    time.sleep(10)
  print 'Done.', task.status()

locations = pd.read_csv('locations.csv')

# Transforms an Image Collection with 1 band per Image into a single Image with items as bands
# Author: Jamie Vleeshouwer

def appendBand(current, previous):
    # Rename the band
    previous=ee.Image(previous)
    current = current.select([0,1,2,3,4,5,6])
    # Append it to the result (Note: only return current item on first element/iteration)
    accum = ee.Algorithms.If(ee.Algorithms.IsEqual(previous,None), current, previous.addBands(ee.Image(current)))
    # Return the accumulation
    return accum

imgcoll = ee.ImageCollection('MODIS/MOD09A1') \
    .filterBounds(ee.Geometry.Rectangle(-106.5, 50,-64, 23))
img=imgcoll.iterate(appendBand)

for loc1, loc2, lat, lon in locations.values:
    fname = '{}_{}'.format(int(loc1), int(loc2))

    offset = 0.11
    scale  = 500
    crs='EPSG:4326'

    region = str([
        [lat - offset, lon + offset],
        [lat + offset, lon + offset],
        [lat + offset, lon - offset],
        [lat - offset, lon - offset]])

    while True:
        try:
            export_oneimage(img,'Data',fname,region,scale,crs)
        except:
            print 'retry'
            time.sleep(10)
            continue
        break


texts: Understanding Mosaic to New Raster treatment of pixel values?

texts:
I'm working with arcGIS.

I have mosaicked 8 rasters with values ranging from -408 to 180. The new mosaic shows pixel values ranging from 0 to 180.Then I compared cell by cell values from the original 8 rasters to the new one, and they match. (The values in the new raster are just integers and decimals in the original rasters, but values in both are very close, rounded to the next integer). However, the negative values are gone in the new raster.

How does this tool treats pixel values from original rasters and target raster?


texts:
Can anyone explain how to install GDAL/OGR with Python on Windows?

I have Windows Vista and I have tried following the information on the website and it does not seem to bind the exe files for me.

Can someone describe the process, including links to the files/folders I will need?

I have now tried to run the gdal setup with minGW, but this has also failed:


texts:
I try to develop some code that loops through a set of geometries (ordered by id)  and deletes from the next neighboring geometry based on the size, but  my syntax seems at the loop seems to be incorrect.

When I run:

WITH firstblocks AS(

SELECT  ST_Buffer(ST_Line_Substring(the_geom, 8*n/length,
  CASE
    WHEN 8*(n+1) < length THEN 8*(n+1)/length
    ELSE 1
  END),4, 'endcap=flat join=mitre mitre_limit=3.0' ) As the_geom

FROM
  (SELECT
  ST_LineMerge(geom) AS the_geom,
  ST_Length(geom) As length
  FROM example_route_single_feature
  ) AS t

CROSS JOIN generate_series(0,CEIL(length/8)::int) AS n
WHERE n*8/length < 1

)

SELECT the_geom FROM firstblocks;

DO
$do$
BEGIN
    FOR i IN the_geom LOOP
        IF ST_Area(i) <= ST_Area(lead(i))
        THEN
            ST_Difference(i, lead(i)) ;
        ELSE
            ST_Difference(lead(i), i) ;
        END IF;
    END LOOP;
END
$do$

I get the error:

SQL error: WITH firstblocks AS(  SELECT  ST_Buffer(ST_Line_Substring(the_geom, 8*n/length,   CASE   WHEN 8*(n+1) < length THEN 8*(n+1)/length   ELSE 1   END),4, 'endcap=flat join=mitre mitre_limit=3.0' ) As the_geom    FROM   (SELECT    ST_LineMerge(geom) AS the_geom,   ST_Length(geom) As length   FROM example_route_single_feature   ) AS t    CROSS JOIN generate_series(0,CEIL(length/8)::int) AS n WHERE n*8/length < 1  )  SELECT the_geom FROM firstblocks;  DO $do$ BEGIN   FOR i IN the_geom LOOP      IF ST_Area(i) <= ST_Area(lead(i))       THEN            ST_Difference(i, lead(i)) ;         ELSE            ST_Difference(lead(i), i) ;         END IF;     END LOOP; END        $do$ returned 0 [ERROR:  syntax error at or near "the_geom"
LINE 26:   FOR i IN the_geom LOOP
                    ^
]

I'm a PostGIS beginner and not so fit with the PostgrSQL-Syntax, so I would be very happy if someone can point whats wrong with the code.

I'm using PostGIS 2.5.1


texts:
Is there a way to add a non-breaking space to CartoDB legends or info windows?

For instance, I have the following legend:

But now I would like to add some spaces to that legend. I used the following approach:

Hello &nbsp StackExchange

Just adding &nbsp &nbsp &nbsp some text...

     ...

As you can see the &nbsp is highlighted in the code:

But it somehow is not translated, and shows up like this:


texts: Performing unweighted overlay on multiple raster layers

texts:
Is there a simple process in ArcMap v10.2 that will allow me to overlay three raster images, whereby any obscured pixels will be deleted in the output file?

I do not want to retain the values of any pixels that are being covered by upper layers, because any overlapping pixels will alter future statistics that I run to determine count.

This is how I envision my final raster:

The bottom layer is a supervised land cover classification raster that has been smoothed over with a 7 x 7 focal statistic (comprised of five classes), which is meant to be my "base" layer.  The middle and top layers are both reclassified images of the bottom layer, and have been smoothed over with a finer focal statistic to retain some detail.

All of these rasters are the same extent and pixel size -- based on Landsat 8 imagery.

Apparently, I am looking for some sort of unweighted overlay tool that simultaneously removes underlying pixels but I have not been able to find this capability in the Toolbox.  I have already tried using the Merge tool, Combine tool, and Mosaic To New Raster tool.

UPDATE

After playing around with the advice given in the comments, I believe using the Con expression as part of an equation in the Raster Calculator would yield the best results.  I would like to build an expression that states "if the value of [middle layer] is 2 then keep [middle layer's] values, but if the value of [middle layer] is not 2 then use [bottom layer's] values."  Then I would perform a similar calculation to stack the top layer over the combined bottom and middle layers.  If there is a way to stack all three layers using one equation, that would be preferred.


texts: Programmatically turn on/off panels with PyQGIS

texts:
In QGIS, with right-click somewhere in the toolbar area panels can be toggled on/off like shown in the screenshot below (the GPS Information Panel for instance).

I need to do this programmatically, so how to toggle on/off panels from PyQGIS? I found that the QgsInterface class has some addDockWidget and removeDockWidget method, but I really couldn't figure out if this is the right way to achieve this ...


texts:
I am working with thousands of tiled (2km x 2km) 4-band NAIP images.  When I diced the images in Erdas, the program corrupted many of the tiles.  The screenshot shows an example of one of the corrupted rasters. You can see there is a small strip of actual pixel values at the top that range from 0 - 255. It is likely that there could be a strip on any of the sides of the image.  The larger black area contains all 0 values.

My attempts to programmatically scan the tiles for a maximum pixel value of 0 or a unique pixel count of 0 failed due to the small area of legitimate pixel values.  This is a simplified approach I have been using:

import arcpy

# input raster data is 8-bit unsigned integer with 4 bands (CIR)
raster = r'D:\temp\4310605_ne_4_2.tif'

p = arcpy.GetRasterProperties_management(raster, "MAXIMUM")

if p == 0:
     print "there is a problem"

What fast and efficient method can I use for checking 4-band tiff files for these corrupted areas?


texts: Using a Join then Expression in ModelBuilder, no data coming up?

texts:
Have a model that converts an excel table to a geodatabase table and converts a shapefile to a layer.  Then joins the geodatabase table and layer together.  After the join I want to create a new layer with specific attributes using an expression.  However when I click on "get values" to build my expression, no data for any of the fields comes up.

Wondering if I'm missing a step after the join?  Or if ModelBuilder doesn't let you select attributes since the join hasn't actually happened yet?


texts: How to calculate the length of each way in a set in Overpass - use of foreach and sum(length())

texts:
I need help understanding the correct use of foreach and sum(length()) in an Overpass query. I've written a query for extracting all relations and ways of a given type within a bounding box, and then, if it finds any, filtering out the ones that don't have ways with highway=footway tags within them (it also extracts parks, but maybe not relevant to the question, apart from explaining why I've named all the sets).

way(8.97266,-79.58028,9.01195,-79.53501)[leisure=park][natural!=wood]->.parks;
(
  rel(8.97266,-79.58028,9.01195,-79.53501)[natural=wood];
  way(8.97266,-79.58028,9.01195,-79.53501)[natural=wood];
)->.tempOther;
if (tempOther.count(ways) != 0)
{
  .tempOther map_to_area->.aOther;
  way(area.aOther)[highway=footway]->.fOther;
  foreach {
    (.fOther; .fOther >;);
    is_in;
    (
      rel(pivot)[natural=wood]->.rOther;
      way(pivot)[natural=wood]->.wOther;
    )->.other;
  }
}
(
  .parks;
  .other;
 );
out geom;

I would further like to exclude relations and ways that only contain short footways, for example, less than 1000 m. I've tried using if (sum(length()) > 1000) {...} (and tried fOther.sum(length())) within the foreach loop, but these fail to do what I think they should do... calculate the length of each individual way (as far as I can tell, it's calculating the total length of all the ways in the set).

Having read the Overpass Language reference and guide, I'm probably still misunderstanding how foreach and/or sum(length()) work.

How do I use them correctly in this situation?


texts: ArcGIS Server extension crashing after some time

texts:
I have developed an ArcGIS server object extension.
It works fine the first hour or so no matter how many times I call it.
However if I do some other work and start calling it after a few hours I get the following errors in the log:

Creating instance of service 'Australia2D.MapServer'.
Failed to initialize server object 'Australia2D': 0x80004005: Unspecified error.
Failed to construct instance of service 'Australia2D.MapServer'. Unspecified error.

I write to the log whenever the soe enters any of my methods and it does not. It crashes somewhere outside my code.

Basically what the soe does it that it serves out subset of rasters from an Enterprise DB. Each time a new raster is requested it is opened and the IRaster is stored in a Dictionary to speed up the next request.

There are no error logs.

I am not using Marshal.ReleaseComObject().
Some of the SOE examples are using this and some are not.

What could be the problem here?


texts: Is it possible to export a GeoJSON with an array field from a joined layer?

texts:
I'm using QGIS to 1) load a shapefile, 2) join the shapefile data with a CSV-file and finally 3) exporting the layer as GeoJSON. Hence, I end up with something like this, where the properties object contains the original CSV data.

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
         "type": "Point",
         "coordinates": [102.0, 0.5]
       },
      "properties": {
         "prop0": "value0",
         "prop1": "value1",
         "prop2": "value2"
      }
    },...

What I really want to do is to store some of the properties in an array, as such:

"properties": {
    "prop0": "value0",
    "arr": [
        "value1",
        "value2"
]

Is it possible to export GeoJSON in such a format using QGIS? If not, does anyone have any other ideas as to how to edit GeoJSON data like that?


texts: Updating Field Mappings Parameter on table selection in Python Toolbox?

texts:
I am building a python toolbox with ArcPy and would like to associate an input feature class parameter with a field mapping. I want the input dialog to work as it does in many built-in geoprocessing tools (e.g. Dissolve), so that once user selects the input feature class the field map populates with the input fields.

I have tried:

fieldMapParameter.addTable(feature_class_parameter.name) in the updateParameters method but throws error ParameterObject: Get attribute: addTable does not exist
Setting fieldMapParameter.parameterDependencies = [param1.name] doesn't seem to have an effect.

def getParameterInfo(self):
    """Define parameter definitions"""
    feature_class_parameter = arcpy.Parameter(displayName="Input_Table", #0
                                    name="fc_settlements",
                                    datatype="DEFeatureClass",
                                    parameterType="Required",
                                    direction="Input")

    fieldMapParameter = arcpy.Parameter(displayName="Primary Settlement Fields", #1
                            name="primary_fieldmap",
                            datatype="GPFieldMapping",
                            parameterType="Optional",
                            enabled=True,
                            direction="Input" )
    fieldMapParameter.parameterDependencies = [feature_class_parameter.name] # Doesn't seem to work
    params = [feature_class_parameter,fieldMapParameter]
    return params

def updateParameters(self, parameters):
    if parameters[0].altered:
       parameters[1].addTable(parameters[0].valueAsText) # Fails with method not found
    return

Full Toolbox

How do I associate these two parameters so they work like field mapping for the native tools such as Merge/Append etc ?


texts:
I would like to do some basic linear referencing (label the distance along a line every 1000 m e.g.) without using plugins or PostGIS or creating points along the line. I tried working with line_interpolate_point (for both rule based labeling and the repetition of single labels) but I don't think this is the correct direction.

Basically, I want to have a simple labeling style which displays geometry-based values every X m.


texts: Incorrect mean calculation with.UpdateCursor ArcPy

texts:
I find this code here, but the results of calculation are not correct.
Field list is correct:
[u'T34TET20160527', u'T34TET20160606', u'T34TET20160616', 'Mean_16']
I calculated mean with Field calculator in the next column in order to check it, but there are not same.

Do you know why?
TableToCalc ='pont'
field_names=[f.name for f in arcpy.ListFields(TableToCalc,'*T34TET2016*')]
field_names.append('Mean_16') #add your mean field to the list

with arcpy.da.UpdateCursor(TableToCalc,field_names) as cursor:
     for row in cursor:
         row[-1]=sum(row[:-2])/len(row[:-2]) #Calculate mean and assign this to the mean field (-1)
         cursor.updateRow(row)


texts: Find all closest points to unique point in PostGIS

texts:
I have two tables valid_collisions (POINT) and intrct_pts (POINT)
I want to create a new table using a nearest-neighbor type of query to find aggregate values from my valid_collisions. In essence I am looking for the geometry of the intrsct_pts and the aggregated values from valid_collisions based off the distance of the compared geometries.
Example of some code I have
SELECT PT.GEOM, PT.CT, PT.ID,
       SUM(VC.PERSONS) AS PEOPLE_INJ,
       SUM(VC.MOTORIST) AS MOTOR_INJ,
       COUNT(*) AS INCIDENT_CT
FROM   PUBLIC.VALID_COLLISIONS AS VC,
       PUBLIC.INTRSCT_PTS AS PT
WHERE PT.GEOM(LIMIT 1)***
ORDER BY ST_DISTANCE(VC.GEOM, PT.GEOM) DESC;

The LIMIT is an error (just there for temp code), but I know that I need to limit the aggregated values and selections by each unique pt.geom. But how would I go about this through using a spatial call like ST_Distance to find all closest vc.geom values to pt.geom records?


texts: Creating empty group layer within existing group layer using arcpy.mapping?

texts:
Referencing this question, Adding new group layer with ArcPy?, I'm trying to use python to add a empty group within a group layer.  So far I'm able to add a empty group within the map document, however I'm getting an error when trying to add another empty group within the existing group layer.  Here is the error message:

  Runtime error  Traceback (most recent call last):   File "",
  line 1, in    File "c:\program files
  (x86)\arcgis\desktop10.1\arcpy\arcpy\utils.py", line 181, in fn_
      return fn(*args, **kw)   File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\mapping.py", line 88, in
  AddLayerToGroup
      assert isinstance(target_group_layer, Layer) and target_group_layer._arc_object.isGroupLayer AssertionError

Here is my code:

        mxd = arcpy.mapping.MapDocument(r"C:\Temp\test.mxd")
        df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
        groupLayer = arcpy.mapping.Layer(r"C:\Temp\Group.lyr")
        arcpy.mapping.AddLayer(df, groupLayer, "BOTTOM")
        targetGroupLayer = arcpy.mapping.ListLayers(mxd, "Group", df)
        addLayer = arcpy.mapping.Layer(r"C:\Temp\Group2.lyr")
        arcpy.mapping.AddLayerToGroup(df, targetGroupLayer, addLayer, "BOTTOM") # error here
        mxd.save()
        del mxd,df

Can you add empty group layers within an existing group (using v10.1 SP1)?  The ESRI documentation (see image below) states that you can add .lyr to groups (maybe not a .lyr group layer?).

Here is the working code (thanks to Jason Scheirer):

    mxd = arcpy.mapping.MapDocument(r"C:\Temp\test.mxd")
    df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
    groupLayer = arcpy.mapping.Layer(r"C:\Temp\Group.lyr")
    arcpy.mapping.AddLayer(df, groupLayer, "BOTTOM")
    targetGroupLayer = arcpy.mapping.ListLayers(mxd, "Group", df)[0]
    addLayer = arcpy.mapping.Layer(r"C:\Temp\Group2.lyr")
    arcpy.mapping.AddLayerToGroup(df, targetGroupLayer, addLayer, "BOTTOM")
    mxd.save()
    del mxd,df


texts: I am working in QGIS 2.1 with the National Land Cover Database 2006 land cover layer. It is a raster layer with pixel values corresponding to different land cover types. I have clipped this raster to a buffer layer around my study site. Now I want to analyze the proportion of each land cover type within the site radius. I can view the frequency of each pixel value in the histogram tab of the layer properties window, but I would like to see the pixel frequency data in a table so that I can perform calculations on it. Is there a way to do this in QGIS?


texts: How to smoothen and dissolve the output image and how to calculate the pixel areas of different classes in Google Earth Engine?

texts:
I am doing supervised classification of 5 different classes, I tried the code suggested in GEE tutorial but I can't able to understand the way they did. I am attaching my code here. Could someone guide me how to calculate the area of different classes? And I want the area of my study region only it is attached as a n asset (tri).

Can we extract our study region (tri - shapefile) from the Landsat image (mosaic)?

Here is the link to my map
https://code.earthengine.google.com/037ea5a6e8240d4d292f0149e1fb996e

Here is my code:

// filter image classification
var image= ee.Image(ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')

  .filterBounds(roi)
  .filterDate('2013-02-02','2017-02-02')
  .sort('CLOUD_COVER')
  .first());
print (image);

var Trivandrum = ee.Image ('LANDSAT/LC08/C01/T1_SR/LC08_144054_20140126');
var Trivandrum1 = ee.Image ('LANDSAT/LC08/C01/T1_SR/LC08_143054_20150122');

var mosaic = ee.ImageCollection.fromImages([Trivandrum,Trivandrum1]).mosaic();
var composite = {
  bands: ['B2','B3','B4'],
  gamma: 1,
  max: 0.80,
  min: 0.04
};

Map.addLayer(tri,{},"triv");
Map.addLayer(mosaic,{min:10,max:200},"mosaic");

var newfc = urban.merge(agriculture).merge(waterbodies).merge(wetland_coastal).merge(forest);
print(newfc, 'newfc');

// Select the bands to be used in training
var bands = ['B2', 'B3', 'B4', 'B5','B6','B7'];

// Sample the input imagery to get a FeatureCollection of training data.
var training = mosaic.select(bands).sampleRegions({
  collection: newfc,
  properties: ['landcover'],
  scale: 30  // should reflect the scale of your imagery, the scale is 30 meters
});
var trained = ee.Classifier.cart().train(training, 'landcover', bands);

// Classify the image with the same bands used for training.
var classified_ = mosaic.select(bands).classify(trained);

var palette = [
  'FF0000', // urban
  '008000', // agriculture
  '0000FF', //  water bodies
  'E3BC1E',   //coastal areas
  '1EE363', // forest

];

Map.setCenter(76.98244, 8.4686, 11);

Map.addLayer (mosaic, composite);
Map.addLayer (classified_, {min: 0, max: 5, palette: palette}, 'Land Use Classification');

var geometry = ee.Geometry.Polygon(
          [[[76.2593, 8.8932],
          [76.4076, 8.0402],
          [77.314, 8.915],
          [77.4897, 8.2686]]]);

Export.image.toDrive({
  image: classified_,
  description: 'Trivandrum_mosaic',
  scale: 30,
  region: geometry,
  crs: 'EPSG:4326',
  folder: 'GEE',
  maxPixels: 10000000000000
});


texts: Keep the same space between lines in OL3 despite scale

texts:
I have 2 lines (in 1 layer - ol.layer.Vector) for 1 road (each line represents different direction). The problem is that when you zoom out then you see only 1 line because they "merge" together as they are so close to each other (cca 20 meters). Is it possible to style the layer in the same way as google traffic layer - so you could see both lines even in small scale?
Zoomed in:

Zoomed out:

This is how it should look like (zoomed out):


texts: Classifying building walls from LiDAR point cloud

texts:
I am currently working on mobile point clouds in order to extract 3D building models in the residential area. I could successfully classify ground, vegetation and building roof lines, but still have an issue with classifying building walls. The walls are not even classified as vegetation, they are still remained as unclassified. I could probably create 3D building from the roof lines by giving height values, but I want to get the entire building points including walls to create multipatch if I could.

I used FME (LAStools) and ArcGIS in classification process as well as viewing and analysis.


texts: What is a correct input format for using qgis.processing.run("gdal:buildvirtualraster") in a python script for QGIS 3.10?

texts:
struggling by creating a virtual raster using a python script in QGIS.  Currently I have a script as follow:

rasters = [r'Z:\Rasters\raster1.tif', r'Z:\Rasters\raster2.tif']
outfile = r'Z:\Rasters\extracted.vrt'
processing.run("gdal:buildvirtualraster",{'INPUT':rasters,'RESOLUTION':1,'SEPARATE':False,'PROJ_DIFFERENCE':True,'ADD_ALPHA':False,'ASSIGN_CRS':None,'RESAMPLING':0,'SRC_NODATA':'','OUTPUT':outfile})

I also have a newline separated rasters.txt file with the contents as:

Z:\Rasters\raster1.tif
Z:\Rasters\raster2.tif

Which I pass to the function as:

rasters = [r'Z:\Rasters\raster1.tif', r'Z:\Rasters\raster2.tif']
outfile = r'Z:\Rasters\extracted.vrt'
rasters = r'Z:\Rasters\rasters.txt'
processing.run("gdal:buildvirtualraster",'INPUT':rasters,'RESOLUTION':1,'SEPARATE':False,'PROJ_DIFFERENCE':True,'ADD_ALPHA':False,'ASSIGN_CRS':None,'RESAMPLING':0,'SRC_NODATA':'','OUTPUT':outfile})

I have also tried passing the contents of rasters.txt to the function as a string in the format:

rasters = r"'Z:\Rasters\raster1.tif','Z:\Rasters\raster2.tif'"

For all of these I am getting "Incorrect parameter value for INPUT" errors.  I also can't seem to find any information on what a correct parameter value for input should be.  I would greatly appreciate any help.


texts: Showing attribute table of selected features using ArcPy?

texts:
Is there a way to show the attribute table of selected features using arcpy?  Maybe in the mapping module?

I would like to be able to select some features and then have an Add-In button that displays the attribute table for just the selected features.

I understand that this seems rudimentary for GIS people, but I'm trying to meet the needs of a specific request for work.

I'm able to do this now using tkinter and only showing the needed attributes, but it's a bit messy and isn't stable.  I know, The stability issue is a known issue with using tkinter in parallel with ArcGIS. So... I'm hoping I can just bring up the attribute table in an automated "neat little button with the company logo on it" way.


texts: Overpass API: Find nearest intersection to lat/lng

texts:
I have a lat/lng pair, and I need to find the nearest street and intersection.
For example: street X, corner with Y.
I have Nominatim and Overpass running.
My approach is:

Get all ways near the lat/lng with Overpass
Calculate distance between all nodes and my lat/lng
Find the two closest to my lat/lng and that belong to different streets
Those should be X, and then Y

Is there any better approach? I feel like I'm wasting CPU and underutilizing Overpass.
Also, I know that X and Y may not intersect at all. This only means they are close, nothing else.
Any help is welcome.


texts: Classification model using NDVI time series?

texts:
I am trying to form a classification model with the help of NDVI time series data. I want to use methods like SVM and Random Forest. I know how to apply these methods to normal data. But, I am confused while applying it to this time series data.
This is a sample of my data:

Here the data have the different type of crops as classes. After that, each column represents the number of days after sowing.
I was thinking of taking these columns i.e. number of days as my random variables and Class as my response variable to apply random forest and SVM classifiers. But I have a doubt if we can do that?? Is it right to apply classification method or there is some other procedure to do it on this type of data?? Please suggest me how can I apply classification method?? Any help would be highly valuable for me.


texts: Latitudes and Longitude Differences Based on Different Datums

texts:
This question refers to the differences in latitude and longitude values between NAD27 and WGS84 for Washington State.
If I create a point using GCS NAD27, and populate the attribute table with lats and longs based on the points datum, and then change the data frame to WGS84 and I populate another two attributes with lats and longs using the data frame's WGS84 spatial reference system, the coordinates are the same.
Why are the values the same?  Shouldn't there be a shift if the X and Y values that reflects differences in the datums?


texts: ArcGIS Server 10 - Query numeric field with string

texts:
In ArcGIS Server 10, I am running queries against several fields in a feature class.  All of them are sent with the format FIELDNAME ='value' (all queries contain the single quotes).

This is failing in cases where have numeric data.  Is there any kind of setting I can use with the map service to allow me to query numeric data with strings?


texts: Plotting points on a map, the json objects retrieved from a database; using WebGL shader code

texts:
I am using an OSM dataset using which I am plotting those points on a map (google maps API). I am plotting these points using webgl shader code (see:  http://psousa.net/demos/webgl/)

The above link by itself is an extension of the Google maps WebGL example :  http://google-maps-utility-library-v3.googlecode.com/svn/trunk/canvaslayer/examples/hello_webgl.html. The difference between the example I quoted and mine is, I have imported the OSM dataset into a postgres database. I am fetching the result (the points that i need to plot) as json objects. After I get the json objects, I am trying to plot these using the webgl shadercode as in the example link above.

The problem is, I am getting the json points (when I do a console.log) but my webgl code does not plot these points. I do not know where I am going wrong.

Here is my code for retrieving the points from the postgres database using jquery, as json objects:

//data.js

 $(document).ready(function(){

     $.ajax({ type: "GET",
                    dataType: "json",
                    url: "/database.php",
                    data: {
                           format: "json"
                               },
                   success: function(data)
                                     {
                      points = data;
                      console.log(data);
                                      }
                   });
  });

The rest are same as in the example link.

Every other function such as resize(), translateMatrix(), scaleMatrix(), init(), loadData(), update() are the same. The only difference is that, while the example plots random points across the map, I am plotting the points I retrieved from my postgres database as json objects. What am I doing wrong here?

Edit : This is the postgres query I used:

SELECT ST_AsText(ST_Transform(way,4326)) FROM planet_osm_line WHERE highway IS NOT NULL AND ref = 'US 101' AND ST_Intersects(way, ST_Transform(ST_MakeEnvelope((-122.302571660907), (37.5601802504715), (-122.436330806712), (37.8074531813721), 94326), 900913));

and the output is in this format:

"LINESTRING(-122.436063557915 37.7998944536137,-122.436172613391 37.7998804703684,-122.437602282165 37.7996967004673,-122.437703791792 37.7996837109204)"

(I am retrieving points within a bounding box)


texts: Minimum requirements (especially screen size) for QGIs - Android

texts:
I would like to try Qgis on Android devices, but before, I would like to know the minimum terminals requirements, and especially the screen size or resolution, what is the minimum screen size do you recommends? At least to don't have any trouble to use it (I tried it on my smartphone... but some menu where larger than the screen size, so I had to stop my smartphone)!
At the same time, if some people have already used it, do not hesitate to give more requirement (memory, cpu, etc.).

Thanks


texts:
I have a map service with an SOE. The SOE works fine and I can call it through http no problem. I would like to call it from an ArcObjects tool but  My code errors when trying to cast the MapServer object to the ServerObjectExtensionManager. I think this is because I get a reference to the map service using http, not a LAN connection. Is there a better way of doing this? I can't use a LAN connection but I have to believe there is still a way to get an SOE.

public IServerObjectExtension InitGeoDataServerFromInternetServer(String url, String serviceName)
    {
        IServerObjectExtension pSOE = null;
        // Create a property set for connection properties.
        IPropertySet propertySet = new PropertySetClass();
        propertySet.SetProperty("URL", url);
        propertySet.SetProperty("user", "user");
        propertySet.SetProperty("password", "pass");
        // Connect to the server and get an enumerator for its objects.
        IAGSServerConnectionFactory2 agsServerConnectionFactory = new AGSServerConnectionFactoryClass();

        try
        {
            IAGSServerConnection agsServerConnection = agsServerConnectionFactory.Open(propertySet, 0);
            IAGSEnumServerObjectName enumServerObjectName = agsServerConnection.ServerObjectNames;
            enumServerObjectName.Reset();

            // Iterate through the objects to locate the map service.
            IAGSServerObjectName serverObjectName = null;

            while ((serverObjectName = enumServerObjectName.Next()) != null)
            {
                if (serverObjectName.Name == "Assessor/Assessor" )
                {
                    IName name = (IName)serverObjectName;
                    //ESRI.ArcGIS.Server.IServerObject pSO = (ESRI.ArcGIS.Server.IServerObject)name.Open();

                    IMapServer pO = (IMapServer)name.Open();

                    ESRI.ArcGIS.Server.IServerObjectExtensionManager pSOEM = (ESRI.ArcGIS.Server.IServerObjectExtensionManager)pO;
                    pSOE = pSOEM.FindExtensionByTypeName("CreateVersionSOE");

                }

            }
        }
        catch (System.Runtime.InteropServices.COMException sex)
        {
            System.Diagnostics.Debug.WriteLine(sex.Message);
        }
        catch (Exception sex)
        {
            System.Diagnostics.Debug.WriteLine(sex.Message);
            //throw;
        }

        return pSOE;
    }


texts: Creating smallest-circle covering polygon using R with sf

texts:
How can I create the smallest-possible circle which covers a polygon geometry using the sf (preferred) or sp packages in R?
One reprex using GADM data:
# Load packages
library(sf)

# Load Nigeria national boundary polygon
nigeria_raw <- readRDS(url("https://biogeo.ucdavis.edu/data/gadm3.6/Rsf/gadm36_NGA_0_sf.rds"))

# Create geometry of the smallest-possible circle
# which covers all points/entire polygon of nigeria_raw

I've not been able to find a solution online.
A solution using sf would be preferred.
Smallest-circle problem: https://en.wikipedia.org/wiki/Smallest-circle_problem


texts: Can not find "application/json;type=geojson" in Geoserver Tile cache format

texts:
I want to server TMS service using geoserver. I know it can be done using GeoWebCache which is loaded by default on geoserver, I have downloaded last version of geoserver and installed it,there must be 'application/json;type=geojson' in Geoserver Tile cache format, but there is not any,for example this tutorial shows it

So, why can't I find it in geoserver admin. Should I do anything else?


texts: QGIS 3.6 - Updating layers from changing excel spreadsheet

texts:
I am building/editing a table and utilising Excel as the vehicle.
When I change the spreadsheet I need to produce a new csv file and then re-import to generate new layers.
Handling Excel spreadsheet in QGIS? gives me some options.

But would it be easier to use a database for the process? I can update/edit the database which could then automatically update the layers? (Postgis, Spatiallite, GeoPackage?)
If so, which database should I be moving to? There will be a steep learning curve for all of them.


texts:
I want to allow a non-GIS user edit a CSV file which is hosted on QGIS Cloud on a PostGIS database. My method is to load the table from the QGIS Cloud database into Excel via ODBC.
Here are my steps so far:

Downloaded the psqlodbc installation file from this link (64 bit).

Opened the control panel | Administrative Tools | ODBC Data Sources (64-bit).
Strangely only 34 bit is available here....

User DSN | Add | PostgreSQL ANSI | Finish (see screenshot)

Input database details as per QGIS Cloud connection FAQ

When testing the connection:


texts: Setting up a Workforce Field - Connecting Workforce and Collector

texts:
I am working on connecting Workforce Developer and Collector and I am trying to set up a specific field as a collector field. There are 3 fields that I can choose from, but none of them are what I am looking for. Are their stipulations to choosing a collector field?

UPDATE:
I was able to get the field that I was interested in. I just had to make sure that the field was a text field. However, when I go to create a work order the work order ID does not populate even though the ID has been set up to correspond with the feature class field. Is there some reason why my workforce ID is not autopopulating?


texts: Problem with geotiff in qgis only in some PCs

texts:
Helllo

I'm having some weird problem opening a 16bits geotiff with qgis. In my pc it opens and displays perfectly, but in another pc (with the same version of qgis) it displays all wrong
the image getting pixels of several colors without any sense.

Any ideas?
Thanks


texts: ArcObjects - accessing layer data table from C#

texts:
I've spent all day yesterday trying to trace all errors but I'm just completely stuck now. Google isn't helpful either.

I'm trying to create a small C# tool which opens a MXD file in background and reads data from one of the datasets (layers). I can happily open the file, I can list all the layers inside with no problems. But no matter how hard I try, my FeatureClass or DataTable (tried both techniques) are always null.

I'm simply trying to get a reference to an ITable or IFeatureClass object to perform a QueryFilter on it.

EDIT: Solution

I really don't know why my approach using IMapDocument didn't work.  Documentation describes it as the best and most efficient way of accessing your data in background. The final solution involved creating an instance of ArcMap, creating a workspace and opening a feature class within it.

    ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Desktop);
    m_application = new ESRI.ArcGIS.ArcMapUI.MxDocumentClass().Parent;

    IWorkspace2 workspace = FileGdbWorkspaceFromPath(@"F:\Projects\!Water\Sewers.gdb");

    IFeatureWorkspace featureWorkspace = (IFeatureWorkspace)workspace;
    IFeatureClass featureClass = featureWorkspace.OpenFeatureClass("Sewers_v5");

And

public IWorkspace2 FileGdbWorkspaceFromPath(String path)
{
    IObjectFactory objFactory = m_application as IObjectFactory;
    Type shpWkspFactType = Type.GetTypeFromProgID("esriDataSourcesGDB.FileGDBWorkspaceFactory");
    string typeClsID = shpWkspFactType.GUID.ToString("B");
    IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)objFactory.Create(typeClsID);
    return workspaceFactory.OpenFromFile(path, 0) as IWorkspace2;
}

It's probably not the prettiest solution (made up from 5 different source codes) but it works.

I'm really disappointed by the ArcGIS documentation. So detailed and rich, and yet so useless.


texts: Openlayers and KML layer projection issue

texts:
I'm not a GIS expert at all, but I've been asked to add a KML layer to an Openlayers map which probably needs some tweaks and that's where I get confused.

Let say I have an Openlayers map with the following projection:

Projection: EPSG:900913
DisplayProjection: EPSG:4326

The KML file that I have contains Geo-coordinates like these:

44.22912497884251,41.93612718554421,0 44.22898876056436,41.93614742383577,0

When I'm trying to show this data on the map, I got nothing, however another KML file shows perfectly fine. It contains some test data like these:

4651136.141979168, 5265323.674468365 4997243.006006268, 5205611.478023042

As far as I can see, there is a projection issue here. What I don't get is the first data looks more correct according to my map projection settings, but it doesn't show up any where or at least I can't find it any where. The second file seems to have different projection, maybe EPSG:3857 (a.k.a. Spherical Mercator), but it will go straight to the map.

We have also ArcGIS installed in one our computers, however as I said, I'm not a GIS expert. If somebody could guide me either through the exporting KML file in ArcGIS (I have the original shape files as well), or the Openlayers settings to show up the KML file, I'd really appreciate that.


texts: Convert Georeferenced imagery to vector data using python

texts: Are there any opensource tools I could use to take a projected georeferenced image such as a jpg with world file, that is UTM projection such as NAD83 and convert it into distinct vector polygons using python or possibly just some software?  The end result would likely need to be a shp file.


texts: How do I trim the first 5 characters of several layernames at once?

texts: Due to a clip-batchprocess and the tool "package layers" I end up with a geopackage with several layernames that start with "clip_". I would like to trim all the layernames by the first 5 characters instead of manually changing all the layernames. How can I do that?


texts:
I am trying to create a DEM from Vector points data extracted from Google Earth. What is the best method?

I get an error thrown each time, while trying to rasterize Vector point data to raster. Pictures attached

What is the guidance for right parameters to be put into the Processing fields of Rasterize (Vector to Raster)?


texts: Geophysical Gravity/Magnetic data CSV to QGIS to Shaded Contour

texts:
I am learning how to use QGIS and I have some geophysical data; aeromagnetic and gravity data that I would like to visualize. The data is currently in CSV format and I am unsure how I should go about processing the data to import it into QGIS. My logic would be to import it as a raster with each point being located with the x and y values, then to interpolate by some means and contour whatever value corresponds to the gravity/magnetic reading for that point before shading based on some colour scheme.
The gravity data for example has column headers with the associated text/data, headers text shown in the screen grab below.

My main questions are:

is my logic correct?
If so, out of the possible options for values to contour with, which one should I use? e.g. observed gravity, bouguer anomaly etc.
Similar question goes for the magnetic data?


texts: How to name a point where we are "back to the initial route"?

texts: I have an initial route, that I leave (for instance in order to refuel). How to name the point where I'm back on the initial route?


texts: Output layers order in QGIS Graphic Modeller (algorithm execution order / dependencies)

texts:
I created a simple model in QGIS graphic modeller that takes in input parameters like drive time ranges and  start location and outputs 3 styled layers including "start_point", "drive_time_zones" and "nielsen_leh_de_2019" competitors. The model functions correctly except that some algorithms need to be run after some other processing algoritms are finished in cases when those are not directly related (one is input to the other one). I tried to handle it with adding extra algorithm dependencies but without success (like adding dependecies as marked in red below in the model schema)

The exact problems are two. Firstly, "nielsen_leh_de_2019" layer comes out on top of "start_point" layer. I want the order to be switch (as marked in blue in the model result map below). Secondly, the final "Export atlas layout to image" algorithm exports a blank OSM map because it does not recognize that it needs to run after all other algorithms are finished.

I tried to go about it using algorith dependencies in many ways but it does not have any influence.
What is the correct way to use algorithm dependencies to have the algoritms run in desired order?


texts: Deleting feature classes from a geodatabase using arcpy.ListFeatureClasses()

texts:
Following up on Using FeatureClassToFeatureClass from list to concatenated list?

The error I get is "Parsing error SyntaxError: invalid syntax (line 10)".

Using the methods from above I am unable to get following to work.

I want to delete selected feature classes using filtered list. This example is a simplified version of what I need.

# Import arcpy module
import os, fnmatch, arcpy
arcpy.env.workspace = r"D:\_StuData\GIS\Projects\PersonalProjects\Product\test.gdb"

# variables:
fcs = arcpy.ListFeatureClasses()
filtered = fnmatch.filter(fcs, '*copy')

# Delete selected fcs
for fc in filtered
    arcpy.Delete_management(fc)


texts: Specify units when buffering arcobjects geometries

texts:
I have a point geometry that I'm trying to create geodesic buffers on.  The input geometries are in a geographic coordinate system and I would like to specify linear units to buffer on.  However I cannot find the correct implementation of the LinearUnitsClass to use.

If I create a buffer object with

IBufferConstruction buffer = new BufferConstructionClass()

I should be able to buffer a geometry input by specifying the distance units. However the input must be an Unit class object and I can't figure out how to implement that.  I'm probably thinking about this the wrong way.

Main objective here is to take a point and find all features (in another layer) that are within 250 feet  -- or to put it another way all features that intersect a 250 ft buffer on the input geometry.  If I get the buffer polygon I can use a SpatialFilter object to pass to my feature search.


texts: Setting properties of ArcGIS Server services over command line?

texts:
I would like to set the pooling property of 100 ArcGIS 10.3.1 map services. This task takes about 2 minutes for each service if I use ArcGIS Server Manager because it restarts the service automatically.

Is there a command line method to do this for all services?


texts: nearest neighbor within table using PostGIS

texts:
using PostGIS 2.0 and PostgreSQL 9.1 I have a table of photographic observations, sco_photos. The location column contains a PostGIS geometry point specifying where the observation was made. The table is not enormous (about 7,000 rows).

I am trying to write a query that will return, for each observation, the nearest other observation in the same table, along with the distance between the two observations. I can easily do this for a fixed point, e.g., how far is each observation from some central point, but the trick is locating the nearest other point in the table and measuring the distance to it.

This is the closest I've come thus far:

SELECT a.photo_id, b.photo_id, st_distance(a.location, b.location)
FROM sco_photos a, sco_photos b
WHERE st_distance(a.location, b.location) =
      (select min(st_distance(a.location, c.location))
       from sco_photos c where c.photo_id <> a.photo_id);

...but this seems just to hang. Performance for this application is not terribly important, as I need this for analysis not for real-time querying.

Thanks for any advice.


texts:
I have a lot of points. In the attribute table I have the coordinates(coordx / coordy).

How can I move all points, so that they are corresponding to the coordinates value in the attribute table?


texts: Is there any similar function to ST_MakeEnvelope from PostGIS in SqlServer?

texts: I need some function which make the same as ST_MakeEnvelope from PostGIS to use on a SqlServer DB. I need to get only some points from a geometry column to make my application run faster.


texts:
I recently found a video on YT, but wondering what map projection the speaker is using while presenting the map at min. 10.54?

I would like to reproduce the results he is showing. Anyone having an idea on how to solve the problem?


texts: GeoServer WMS removes color from 4 band image

texts:
I have loaded a set of 4 band (R,G,B,NIR) GeoTIFF files in GeoServer and I am accessing the data via WMS. My problem is that whenever the NIR band contains a value of 0, that in the WMS output image the R, G and B bands are also cleared of their values and set to white. This is not what I want.

The input GeoTIFF images have no nodata value set. So why is GeoServer use the NIR 0 value as nodata? Below are the settings:

I tried to use the input and output transparent color, but once I set them the WMS only returns 3 bands in the image and the 4th NIR band is no longer there, so that doesn't help either.

And how can I configure GeoServer to just return the data in the input image without changing it?


texts: Inserting union of two tables into another table using PostgreSQL?

texts:
I am working on a research project. I have two tables, one of a polygon that represents a severe weather risk area, one that represents the US. I want to merge them into one output polygon. I am stuck with the syntax of how to do this.

spc is the table of the risk area, simple_us is the us. I want to output the union of these two table into the geometry column of a new table called spc_clip.

My current command is:

INSERT INTO spc_clip (geom) VALUES (ST_Union(a.geom, b.geom)), spc as a, simple_us as b;

As it sits now I get a syntax error at spc when I run the command.

Any ideas on how to proceed?

Next I want to export the geometry of the union to a JSON. I am not sure how to do this. I cannot test this line until I get the above union done.

SELECT row_to_json(t) FROM (SELECT name,geom FROM spc_clip) t;

NOTE: name is already in the table and works fine. I am just struggling with the join and the json export.


texts: Sum area of specific poly attributes within raster grid cell

texts:
I have a polygon layer (A) containing several different land use types (road, building, etc).  I also have a different GIS layer (B), a raster of NLCD land use/land cover.  I would like to sum/calculate the area of certain attributes in A (say roads) within each grid cell of B, to use in a comparison. Both layers are in the same projection.

How does one go about doing this in ArcMap?

Edit: Here are pics to help understand:

The polygon layer:

The polygon attribute table:

The raster layer sitting on top of the area of interest:

The end result I desire is another raster, with the same footprint as my red/pink raster, that has a value of the area of the 'road' polygon within each grid cell.


texts: In TileMill/CartoCSS how do show several markers for the same point?

texts:
I'm creating a custom map using TileMill to show local walks.

I want to show Points Of Interest along the path.
I've exported the points from Google Maps using KML and I want to show an icon on a background for each.

I have a pin icon and the icon for each location but I can't work out how to display one on top of the other.

I've looked into using a point and a marker but point doesn't allow color styling and I would like to color the background pin as well as the icon.


texts: ArcMap line symbol with two different ends

texts:
Is it possible to create a line symbol that has a round line end at the one end and a straight line end at the other using ArcMap?

I am only able to do something like this with arrow line symbols:

<----

---->


texts:
I want to use the font Air Soft Condensed heavy in Tilemill. I placed the font file Positype-AirSoftCond-Heavy.otf into the font directory of my project:

/home/<user>/Documents/MapBox/<project>/fonts/

The palette.mms file has the line

 Map { font-directory: url(./fonts); }

How can I reference now my font to use it?


texts: OpenGeo/GeoServer WFS - getting only as much data as is needed

texts:
There's probably a term for this that I'm missing that would help me find the answer...

My question is similar to a previous question of mine. I am using GeoServer (as part of the OpenGeo stack) as a WFS server. Currently when I load my map it downloads the entire layer - roughly 10MB per layer, roughly 10 layers...it takes a while.

I would like it to work in the way that WMS does - only download enough to show what is needed.

I need to use WFS rather than WMS as I need to be able to click on polygons, edit their attributes etc.

Is there a word for this?

Is it possible with GeoServer/OpenGeo/GeoExt?


texts: How to transpose and summarise attributes in FME?

texts:
I have a table information and that needs to be transposed like shown below screen shot.

Parcel number is maybe more then 1000 and each parcel bldg type more then 30 types.

Please can advise to achieve this task.


texts:
I am trying to calculate summary statistics for a raster file. I am encountering a problem when loading my raster into Python.

When I load the raster and I plot it, it is flipped 180 degrees. This means it is upside down.

When the raster is loaded into QGIS or ArcMap is properly displayed. How can I proceed?

import os
import sys
import glob
import rasterio
import gdal
from rasterio.plot import show
from rasterio.transform import Affine
import networkx as nx
import geopandas as gpd
import osmnx as ox
import matplotlib.pyplot as plt
from rasterstats import zonal_stats

os.getcwd()

path = "path_to_data"
gis = "path_to_vector"

raster = "path_to_raster**"

raster_yearly_path = [f for f in glob.glob(raster + '*ppty.tif', recursive=True)]

raster_yearly = []

for raster in raster_yearly_path:
    local = rasterio.open(raster)
    raster_yearly.append(local)

ax = show(raster_yearly[0])

I cannot move forward since the basins will not be in the proper extent. Please check the image below. I need to process more than 500 rasters. This is why I am using Python.


texts:
When I am using the raster package to predict a surface probability map of species distribution using the following functions in R, the following error occurrs:

CA132122200=raster("F:/CBG maps-UnProj/anlytic final map/mov472/landmergbackg1.tif_mw2/ca_132.tif")
mod=gam(y~s(CA132122200) ,data=training, family=binomial, gamma=1.4, logLink="logit")
 prediction <- predict(object = CA132122200, model = mod, type='response', progress='text')

errors:
object is not a matrix:
In addition: Warning messages:

1: In predict.gam(model, blockvals, ...) :
  not all required variables have been supplied in  newdata!

2: 'newdata' had 2032869 rows but variables found have 2315 rows"


texts: Separating parcel polygons that share boundaries in ArcMap?

texts:
I have parcels in one layer which are simple polygons. As one parcel is against it's neighbor the lines of that boundary are shared.  This makes it very difficult to symbolize the parcels in different outline colors.

I found a way to deal with it here: Displaying adjacent polygon lines beside each other instead of overlapping

Unfortunately, this fix does not hold when I export into KMZ using the Map to KML tool.

Can someone help me figure out a solution for this that can give me a good solution in the end product of a Google Earth KMZ file?

This is what I am actually dealing with in GE.


texts: Performing buffering and intersection with script in PyQGIS

texts:
I wrote a script that performs the intersection of two layers:

The first Layer -> layer1, represents the rivers
The second Layer -> layer2, represents a set of circular dots distributed in a geographical area (at each point corresponds to a numeric value that exceeds a threshold value)

The result of the intersection should be highlighted with an operation of buffering.

The new geometry obtained with Intersection + Buffering should produce a layer that surrounds the traits of the rivers and the point or more points that intersect.

Writing from console PyQGIS I have not encountered any errors. To see the new geometry I added new layers to the previous ones (from "add vector layer") but I have not seen anything new.

This is the code of my script:

from qgis.core import *
from qgis.gui import *

canvas = qgis.utils.iface.mapCanvas()

layers = canvas.layers()
layer1 = canvas.layer(0)
layer2 = canvas.layer(3)

from qgis.analysis import *

overlayAnalyzer = QgsOverlayAnalyzer()
overlayAnalyzer.intersection(layer1, layer2, "/home/..../Scrivania/..../output.shp")

layer3 = canvas.layer(0)

geometryanalyzer = QgsGeometryAnalyzer()
geometryanalyzer.buffer(layer3, "/home/..../Scrivania/..../buffer.shp", 500, False, False, -1)


texts: Selecting lowest values in from attribute table in ModelBuilder?

texts:
I am trying to select the lowest value from an attribute table within model builder. I have used the Sort tool to sort the values in ascending order but when I add the Select tool and go to enter an expression there is no value selectable. From my understanding, the Sort tool assigns an FID to each value with the lowest value being 1 in which I could use the Select tool and enter the expression: "FID = 1" to get my desired output.

I am relatively new to ModelBuilder.

Any suggestions on how to work through this?

Current Model


texts: QGIS how to save loaded min/max values for raster layers?

texts: I am working with many layers.  Each time I open my saved project, I must reload the min/max (properties/style) values.  Is there a way to turn off the default 2% cut?


texts: Opening up a QGIS project file on a different computer

texts: I've tried doing this but when I save the QGIS Project to my desktop its only a few KB's so when I try to open this up on a different computer all the information inside the actual file is missing.  I'm wondering if there is a way to save all the data to a location or how to find where the data is currently being saved.


texts: Correcting effect of solar elevation angle in tagged areas for Landsat 8 L1T images using ENVI

texts:
Someone said the L1T images did not correct the effect of solar elevation angle.
I want to know how can I get the parameters for doing such correction in ENVI?


texts: How to convert drawn vector feature from EPSG:900913 to EPSG:4326

texts: My script allows you to drawn shapes on an OpenLayers map, and then to get the WKT information from that shape, except the returned coordinates are in EPSG::900913 and I need them to be in lat/lon.


texts: "column does not exist" for SELECT ST_Point

texts:
I am using psql (9.5.2) and postgis 2.2.2. In the psql shell,
I have created a spatial database and connected to it.  I want to insert a point into the postgis geodatabase.

To do so, I first create a table in the database:

CREATE TABLE table1 (
  p_id INTEGER PRIMARY KEY

);

Then I want to create a point using:

select st_point(101, 20);

However, I got an error:

ERROR:  column "101" does not exist

What is the reason of this error?


texts: Why does the composite raster pixel values not cover the full range of RGB and what does this mean for reclassification after conversion to int?

texts:
I have about 30 or so composite rasters in IMG format that I need to 1.) Preserve the date that the raster was created and 2.) Combine all of them. My workflow  is:

Step 1: raster calculator, multiply floating raster by 10 so that decimal values are not lost in integer conversion
Step 2: convert floating raster to integer
Step 3: create attribute tables for the rasters
Step 4: Add a field
Step 5: Calculate field to add the date which will be taken from the name (the datetimes are in the names)
Step 6 create raster mosaic - make sure it is still an integer.

This seems easy enough, but I'm not the best informed when it come to rasters. After step 2 (or step 1 if I multiply it by an integer) I get the raster with an attribute table, but it has lost its RGB coloring scheme. I thought that perhaps I could assign RGB based on pixel values as a step 7 once everything is done (somehow) but the pixel values do not quite match the RGBN individual bands. The range of the values of the new integer dataset is 3760-16180 (meaning the range was 376-1618 before Calculate Raster), which are the same values as the first band only (Red I think) and the composite stretched value range. Why is this? Why does the composite band  use only these values as well when there are values lower and higher than this in the other bands? How can I classify the pixels as RGB if the full range for each is not there? Also, if you have alternatives or suggestions to my workflow, please let me know.


texts:
I have 1000+ irregular shaped polygons where I would like to determine the orientation (azimuth) of their long axis. The polygons overlap in a few cases, if that matters. I am using the Minimum Bounding Geometry (Data Management) tool with the "Convex Hull" geometry type. When I run the tool, the resulting attribute table does not list anything except the length, area, and original feature ID -- not the MBG_Apod information and orientation  as described here.

Why is this occurring and how can I fix it?


texts:
I am trying to decode a .FIT file from a Garmin Instinct GPS Watch. Sourced from the file Lctns.fit.
However, I am unsure how to interpret the entries values as coordinates.
Given this entry:
# Output from parsing the FIT file.

* unknown_1: 571162134
* unknown_2: 92300401
* unknown_254: 12
* unknown_3: 63
* unknown_4: 4930

# The same entry as Lat / Long

N 47.52.4554
E 007.44.1918
Elevation: 486m

How can I determine what coordinate system is being used inside the FIT file as to convert to Latitude and Longitude?
I have tried converting from UTM with no success: http://www.rcn.montana.edu/Resources/Converter.aspx


texts: Using CAD data that is not to be projected in ArcMap?

texts:
I need to incorporate CAD data into my geodatabase but the usual steps seem to require defining a projection or projecting the data in some way.

All of my geodatabase's feature classes are in GCS_WGS_1984 which is technically not projected, so I want the CAD data to be in GCS_WGS_1984 as well.

Is there a way to georeference the CAD data and convert (ie CAD to GIS) it in this case?


texts: Relations and working with tables in QGIS

texts:
I am working on a project where I needed to study soil types at my research location. I am using SOTER. Just for record, I am working in QGIS 2.6 Brighton.

What I have ?

A shapefile (with 25 features). Each of these features have a linkfield column. this column has id that can be used as foreign key to SOTER database
SOTER database with multiple tables: Terrain, TerrainComponent, SoilComponent and Profile are the ones I  am interested  in.

Data Structure::

Terrain table has 17 fields with SUID 1 to 17. "linkfield" of shapefile will be linked to this

Terrain Component table has more than 50 fields, each SUID can have 2-4 terrain components - represented by TCID column

Soilcomponent table assigns at least 2 to 4 different components to each TCID - represented by SCID. this table also has PRID in each row. Thus, each SCID of each TCID of each SUID has 1 distinct PRID that can be linked to Profile table

Profiles table has WRBC column that I am intersted in.

What have I done?

Created relations in QGIS with following parameters

linked SUID of Terrain to linkfield of shapefile
linked SUID of TerrainComponent to linkfield of shapefile
linked SUID of SoilComponent to linkfield of shapefile
linked ProfileID (written as PRID) of Profiles table to PRID of SoilComponent

Now when I click over any feature in shapefile, (with auto-open form) I can navigate to find a column called: "WRBC" of Profile table which has the actual data that I need. I have attached a sample form of a feature of shapefile.

What do I need ?

Now I need to be able to extract all the types of PRIDs that are referred by the "linkfield" from shapefile. And I need to save this into a new table.

The format of the table could be like:

linkfield    PRID     WRBC
1            NP111    its value
1            NP124    its value
1            NP214    its value
2            NP623    its value
2            NP155    its value
2            NP124    its value
3            NP214    its value
3            NP623    its value

and so on...

I need this table for manual study and then later I need to reclassify this data on the basis of the value in WRBC column.

Am I doing anything wrong?
If not how can I proceed?


texts: I am newbie to osm/nominatim. I set up everything correctly and loaded Florida pbf data and geocode /searches were working fine. Today I added data for Mississipi using the setup.php script with --import option and looks like the data got loaded with no errors. However I am not able to geocode any of the mississipi addresses. The old Florida addresses still geocode fine. What am I missing? What should I do after loading new osm states data? I must be missing a simple step. I am happy to post additional information as needed.


texts: How to change QGIS install paths on Windows?

texts:
I will be installing QGIS on client machines soon, and we need to have the .qgis2 folder on the D:\ and not to C. (Windows 7)

Is there a way at install time to have it moved to a specific location?

Or, is there a way to delete the installed locaiton and change the path on startup?


texts: How to Normalize or rescale the NDVI raster image

texts:
I have processed my landsat band 3 and band 4 of different years to derive NDVI values the possible values of Normalized difference vegetation index is -1.0 to +1.0. The bands were atmospherically corrected by using dos1 correction.

NDVI was calculated for my study area and the values that have shown in statistics parameter of every raster layer property as minimum, maximum, mean and standard deviation respectively for each corresponding years are as follows:

1993    -0.15   0.73    0.22    0.099
1998    -0.41   0.85    0.33    0.13
2002    0.0009  0.629   0.165   0.07
2009    -0.061  0.55    0.19    0.068
2011    -0.067  0.596   0.165   0.064

My question here is that how can I normalize all the rasters of each year so that i could see the relative change in each year. In other words, I want to derive the relative index across the study area but not clear that which range should I chose for my case into which my raster values clustered equally interval index. How can I rescale these values

Next,Green vegetation fraction (GVF) was to calculated by using the same equation as Hu and Jia, 2010 used for city of south China , Guangzhou as mentioned equation-8 in paper also named as linear mixing model.

  GVF = NDVI - NDVImin / NDVImax - NDVImin
  where NDVImin and NDVImax has been chosen as 0.05 and 0.70 respectively and these
  values that were derived by AVHRR 10-day composite data.

I want to normalize (or say transform if i am not wrong) all the values to 0-1 so as to compare the fraction of each pixel in green vegetation to examine actual role of vegetation. How can I derive the range for my region.
Do I need to rescale it 0-255 first? Because this is a discrete data not a continuous ranges from -1 to +1.

Linear Spectral Unmixing by using the image based method of collecting endmembers were used in different study to derive the range of NDVImin (bare land) and NDVImax(dense vegetation).


texts: Using ArcPy, I can't panToExtent to intended area using Decimal Degree coordinate

texts:
I'm trying to programmatically pan to a coordinate before I export a map and I've started off by using Identify to get a test coordinate in Decimal Degrees:

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
x = 30.317459
y = -97.694778

extent = arcpy.Extent(x, y, x, y)
df.panToExtent(extent)

arcpy.RefreshActiveView()

When it pans I then use Identify again to check the coordinate and it's in a completely different region. I only have one layer in this mxd and it is in UTM Zone 14 along with the Data Frame.


texts: Geography for Sql Server STIsValid() getting error

texts:
I am getting the following error
Could not find method 'STIsValid' for 'Microsoft.SqlServer.Types.SqlGeography' in assembly 'Microsoft.SqlServer.Types' on the following query.

select ogr_geog from dbo.pod_1730_essential_energy
where ogr_geog.STIsValid() = 0 and ogr_geog is not null

All my other Spatial methods work. The ogr_geog column is a geography type.
Does anybody know the root of this error?


texts: pgr_dijkstra not returning anything when using the same vertex

texts:
I'm using pgRouting 2.3 and in the description of the pgr_dijkstra function states that:

  When the starting vertex and ending vertex are the same, there is no
  path. The agg_cost the non included values (v, v) is 0

My problem is that I need to use pgr_dijkstra One to many

pgr_dijkstra(edges_sql, start_vid,  end_vids, directed:=true)

And sometimes, the same start_vid node could be in end_vids. Yet pgr_dijkstra returns nothing for that pair of same nodes.
Do I have to validate if start_vid is inside end_vids before querying pgr_dijkstra?


texts: Geopandas equivalent to select by location

texts:
I have 2 geodataframes:

A single multi-part polygon
Many points

I have read them in with the following:

poly=gpd.read_file('C:\\Users\\srcha\\Desktop\\folder\\poly.shp')
points=gpd.read_file('c:\\Users\\srcha\\Desktop\\folder\\points.shp')

I want to run the equivalent of 'select by location' in ArcGIS with geopandas and shapely tools to select points that are within the polygon.

How should I go about this?


texts:
I have generated Kernel Density Rasters from GPS species locations. I want to generate percentage isopleth polygons to enable me to identify core areas. I am using ArcMap10.6 and therefore cannot use GME. I'm working on my Master's project and would prefer to use Arc as I am not confident in R.
Is there a simple method to turn my Kernal Density rasters into percentage isopleth polygons?


texts: QGIS 3.x create polygon in specific size - QAD plugin not compatible

texts:
I've noted at Polygon in specific size that the QAD plugin can be used to create specifically sized polygons, however it is not compatible with QGIS versions newer than 2.9.

Any suggestions for a workaround?

I'm presently incrementally resizing my polygons then using the field calculator to approximate their sizes: however, since there are many polygons and they have to be very specific sizes, this method is prohibitively time consuming.


texts: Comparing different degrees of local spatial association

texts:
I have two data sets and want to compare/evaluate their degree of spatial association expressed in terms of the LISA statistic.

The local Moran's I statistic gives for each point a (pseudo) p-Value and a correlation coefficient. I want to make a statement like "The points of data set a are more locally correlated than the points of data set b". To do so I though of simply comparing the sum of the absolute values of the significant (p<0.05) correlation coefficients.

Do you think this is a valid approach? Alternatively it might be also reasonable to form the squared sum in order to give more weight to large coefficients.


texts: QGIS Lat/Long Postcodes Coordinates Being Placed Out At Sea?

texts:
I've read plenty of posts on this site regarding coordinate systems but none of them have helped solve this issue. I'm using Ordinance Survey Lat/Long coordinates for postcodes but a small number are being placed out at sea above Northern Wales (see image).

Surely if the coordinates were completely wrong, it would be more than a fraction of the total points that were incorrectly geolocated.

All layers are in OSGB 1936 / British National Grid (EPSG: 27700). Changing this to WGS 84 doesn't resolve the issue.


texts:
I'm attempting to upload a KML file to Google My Maps.  But it says that the file I'm attempting to upload is either non-parse-able, or too large.

The KML file in question is generated by a program I've written.  I've tried removing all of the polygons save three and the KML parses and uploads correctly.  So I'm fairly sure it's not an error in my KML syntax.  I've also validated it in

Here is the output from my program:

POLYGONS WRITTEN: 53
TOTAL POINTS WRITTEN: 4418
SELECTED OUTPUT FILE: output.kml

I run:

  [foo@localhost test]$ ls -lh

  -rw-rw-r--. 1 foo foo 108K Apr 22 16:48 output.kml

and can see that the file size of output.kml is 108 Kilobytes.

According to Google, the max file size is 10MB.   So I'm well within the limit - memory-wise.

I don't understand why my file is being rejected.  When I remove the majority of the polygons, the file is accepted.  So the syntax is fine. But I seem to be within the file size limitations.

What's going on here?  And how do I get around it?


texts: QGIS, Query out some attribute columns temporarily

texts:
I am looking for a work-around to remove some fields (attribute columns) from a layer's table by a query.

For example, there are 10 columns in a table, remove 9 and leave only one (virtually).

Now, if I run 'Join Attribute by Location' tool, I should only get one column attached to the new 'join layer'. This is so that I could get only the desired attribute column in the output 'join layer'.

Once I finish this task, I should be able to display all original columns back as they were in the layer.

I know there is an equivalent in ESRI's ArcGIS, but not sure about QGIS.


texts: Token based authentification of QGIS Server services

texts:
In ArcGIS for Server a token for a service can be generated in order to restrict the consumption of such service to one single url (Such approach might be necessary due to some strange legal framework in Germany e.g.).

Question: Is there something equivalent available in QGIS Server?