texts: Intersection between polygons layers is empty

texts:
I'm trying to compute the area of green surface for cities with QGis from open data found online. I have polygons on a layer for the cities boundaries and polygons on another layer for the parks.

I first decided to clip the polygons of the parks but as the boundaries of the cities are adjacent, the parks are not cut along the boundaries.

I then chose to intersect both layers as it would cut the parks alongside the boundaries. But the result is empty! Any clue?

I've already checked the CSR which are the same for both layers.

For example, in this picture, there are two parks (North one foret communale du canton du Cerisier and South one foret communale) and two cities (376 West and 519 East). As the limits of the parks are not the same as the ones of the cities, it should return 4 green areas: two for the North forest (one for each city) and two for the South forest (one for each city).

Here are the files I'm using: https://files.fm/f/cpuyqu2f.

Edit: Seems like the problem comes from park layer. I tried with new layers and the polygons intersection works.


texts: Compute convex hull of points in a feature collection

texts:
I uploaded an asset (which was originally a shapefile of points) and imported it into my script as a feature collection. I would like to compute the convex hull of the points, but it seems this function is only available for geometry objects. I have been trying several methods to convert the feature collection to multi-point geometry, but nothing seems to be the correct incantation. How would one convert a feature collection to multi-point geometry (or am I missing the mark in how I'm approaching this)?

I'm open to any GEE solution that accomplishes the main goal.


texts: How many bookmarks can you have in QGIS?

texts:
I am currently working proactively for a big upcoming project at work, where I am utilizing the bookmarking function in QGIS. My problem is that after a certain extent my "save new bookmark" function does not create a new bookmark for the specified location; instead it marks an already existing bookmark and wants to overwrite the name of said bookmark.

Currently I have 256 bookmarks in my project and I will need at least 4-5000 more in order to have all the specified locations needed for the project.

Is there a way around this "invisible cap" or is the only solution to export my 256 bookmarks, delete the list from the project, create a second wave of bookmarks (rinse and repeat) and import them as needed - or is there a more elegant solution to my problem?


texts: Access icons from within QGIS plugin with relative paths?

texts:
I have constructed a QGIS plugin using QGIS Plugin Builder 3.0.3.

Within my plugin, I need to load a .png file as a Qicon object.  What needs to be done for the plugin script to locate the .png file when it is run within QGIS?

For my situation, my .png file is called line_1x.png.

I have modified the resources.qrc within the plugin to show a reference to the icon as shown below:

<RCC>
<qresource prefix="/plugins/geomAttribute" >
    <file>attribute_table.png</file>
    <file>line_1x.png</file>
</qresource>

I compiled the resources file using the following command pyrcc5 -o resources.py resources.qrc.

Within my plugin python script I try to construct the QIcon object using:
icon = QIcon('./line_1x.png')

However, the plugin python script can't access the line_1x.png file.  When I use a full file path to the icon file instead of a relative path the icon image file is read properly.  Hence, using the following does work:
icon = QIcon('C:/Users/pippi/AppData/Roaming/QGIS/QGIS3/profiles/default/python/plugins/geomattribute/line_1x.png'

How do a I get a QGIS plugin to read an image for an icon using a relative path?


texts:
I am trying to make a simple map with the following legend:

Whatever I change in the settings I can't get the other boxes of the predicted deprivation section to be place under the group title. Ideally I would like to have 4 columns UNDER the group title with 2-2-2-1 boxes. Any ideas?


texts: Why general.runalg("qgis:linestopolygons"...) work only over first feature?

texts:
I have a standalone script that call general.runalg(qgis:linestopolygons...)  all ok, but only convert the first feature, and I do previous "selectAll", why ?

Full code:

#!/usr/bin/env python

import qgis
from qgis.core import *
import sys, os

app = QgsApplication([],True, None)
app.setPrefixPath("/usr", True)
app.initQgis()
sys.path.append('/usr/share/qgis/python/plugins')
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *
Processing.initialize()

inFile="temp.shp"
poly_output_file=os.path.splitext(os.path.basename(inFile))[0]+"_poly.shp"
vlayer = QgsVectorLayer(inFile, "temp", "ogr")

print "Elements loaded:"+str(vlayer.featureCount()) # shows number elements loaded

vlayer.selectAll()
general.runalg("qgis:linestopolygons",vlayer,poly_output_file)
# note
# you can also call (but same result), only 1 feature into output:
# general.runalg("qgis:linestopolygons",inFile,poly_output_file)

# Exit applications
#app.exitQgis() # give segmentation fault
app.exit()

The temp_poly.shp file is created, but only have one polygon (the first feature).

If I try to convert the same file in, but using QGIS desktop GUI all features are convert to polygons, so is not a problem from source file.


texts: Batch Convert Folder of GPX to SHP in QGIS

texts: I have a few GPX files (up to 35 Mb, with a few tracks each one) that I want to convert to SHP. QGis is faster and easiest, than ArcGIS to convert GPX tracks to SHP, but is there a way to batch process it and keep the same file name for each new SHP?


texts: How to interpret a WKB (Well-known-binary)?

texts:
I read this example of a WKB (well know binary) from "PostGIS in Action" 2ed.:

SELECT ST_AsBinary(
ST_GeomFromWKB(
  E'\\001\\001\\000\\000\\000\\321\\256B\\312O\\304Q\\300\\347\\030\\220\\275\\336%E@',
  4326
  )
);

I understand the E'..' is an escape string, and the escaped backslash is a separator. But I don't understand the rest of the well-known binary.

Can anyone explain what are the numbers and letters?

They do not look like hexidecimal according to wikipedia consider e.g. 304Q. Also the 336%E@ looks weird.


texts:
I have a shp and I would like the exploded view with geopandas in several polygons,I cannot use ARCGIS and keep that the polygons which have a given.
import csv
import glob
import pandas as pd
import numpy as np
import time
import datetime
import sys
import geopandas as gpd
from shapely.geometry import Polygon

gdf = gpd.read_file("dossier/essais.shp")

time.sleep(10)
xmin,ymin,xmax,ymax = gdf.total_bounds

lenght = 2065
wide = 3069

cols = list(range(int(np.floor(xmin)), int(np.ceil(xmax)), wide))
rows = list(range(int(np.floor(ymin)), int(np.ceil(ymax)), lenght))
rows.reverse()

polygons = []
for x in cols:
    for y in rows:
        polygons.append( Polygon([(x,y), (x+wide, y), (x+wide, y+lenght), (x, y+lenght)]) )

grid = gpd.GeoDataFrame({'geometry':polygons}, crs = gdf.crs.to_string())
grid.to_file("dossier/grid.shp")

#df4.head(1000).to_file("dossier/"+file)
gdf.to_file("dossier/grille.shp")


texts:
I successfully ran a function over a smaller test raster stack (48 layers) in about 30 minutes but cannot run it over my large stack (1462 layers). The last one I tried over the large stack ran from 9PM to 7AM and had not finished processing. I stopped it and tried a simpler function and I have not been able to get it to finish yet.

Large stack information: layers = 1462, total elements = 1275602310, 16.4 mb, ncols per layer = 1405, nrows per layer = 621.

# Import libraries ----
library(tidyverse)
library(dplyr)
library(magrittr)
library(raster)

# Call the source file for the function
source("analysis_functions.R")

# Identify file path for files to stack
t_dsn <- "filepath_here"
p_dsn <- "filepath_here"

# List files to stack
t_list <- list.files(t_dsn, pattern = "tif$", full.names = TRUE)
p_list <- list.files(t_dsn, pattern = "tif$", full.names = TRUE)

# Create stacks
t<- stack(t_list)
p<- stack(t_list)

# Calculate a pixel-wise mean value through the stack
tm <-calc(t, mean)
pm <-calc(p, mean)

# Calculate a pixel-wise Nash Sutcliffe Efficiency value through the stack
nse_t <- myNSE(sim = t, obs = p, obs_mean = mp)

#Lastly, here is the function as defined in the source file
myNSE <- function(sim, obs, obs_mean){
  nse <- 1 - ((sum((obs-sim)^2)) / (sum((obs-obs_mean)^2)))
}

The test stack and the large stack are comprised of the same exact files, just a smaller set.

The mean calculation for the test stack takes about 2 minutes.

The mean calculation for the large stack takes about 30 minutes.

The NSE calculation for the test stack takes about 2 minutes.

The NSE calculation in the large stack hangs for more than 9 hours and I have not been able to complete.

Lastly, I have successfully run all commands until I get to the NSE one. This is the one that seems to lag.

Any suggestions?


texts: raster::focal fails on a selection of invalid row in R

texts:
I try to compute the proportion of 0 in an area of 30km around a given pixel.
The raster has a resolution of 100m and it is the map of France.

For this, I use the function focal() of the raster package (package version: 2.9.23). The analyse works well on a subset of the raster but not on the full raster. The error message is the following:

> Error in .local(x, row, nrows, ...) : validRow(x, row) is not TRUE

Looking for the .local function:

getMethod("focal", signature = "RasterLayer")

I find .local() but did not find its source on
github or in the source code.

validRow function

The function validRow() checks that a number of rows is valid regarding the
raster.

r <- raster()
validRow(r, c(-1, 1, 100, 10000))
# -1 is obviously not valid and 10000 neither because the nb of rows of the
#raster is 180

I did not succeed to identify what operation in the focal() function can lead
to this error.

The problem could be related to this open issue in github.

 Trying to reproduce the error

The most annoying part of the problem is that I do not succeed to reproduce the
error with a minimal reproducible example:

library(raster)

# Get France raster at 100m resolution
p <- getData(name = "GADM", country = "france", level = 0)
p <- spTransform(p, CRSobj = CRS("+init=epsg:2154"))
r_total <- raster(p, resolution = 100)

# Crop and put dumb data
r <- crop(r_total, extent(98900, 150000, 6100000, 6150000))
nbcell <- ncell(r)
# Set values 0 or 1:
r <- setValues(r, rbinom(n = nbcell, size = 1, prob = .9))
# Put some NAs:
values(r)[sample(1:nbcell, size = .9*nbcell)] <- NA

# Compute the focal
fw <- focalWeight(r, 30000,"circle")
test <- focal(r,w = fw, fun = "sum", na.rm = TRUE, pad = TRUE, padValue = NA)

If I could reproduce the error, I hope that It will help to identify the problem
with the real dataset.


texts: ArcGIS Online and Nearmap Basemap Imagery

texts: I am trying to get an exported Nearmap imagery as my basemap for an ArcGIS Online map. I have the Nearmap Image and it's respective georeference file (JPG File and JGW File). Is there a way to upload this to ArcGIS Online and use it as my basemap?


texts: Create polygons of the boundary of a Raster in R

texts:
I want to create two polygons.

One is of the rectangular extents of my raster. I know I can obtain the extent of the raster using
r <- raster("band5.tif")
e <- extent(r)
plot(e) gives me this:

How can I create ONE polygon of the boundary of the raster as shown below?

This is what my raster looks like:


texts: How to select one field and delete all others QGIS graphical modeler

texts:
In the Graphical Modeler in QGIS I have the algorithm to Drop Fields, which you select the fields you want to drop:

However, the vector layer I'm using has many fields but I just need one in the end.
I'm wondering if there's any way to trick that in order for the user to just select this one field as an Input and all the others will be deleted automatically.

EDIT:
Ass suggested by @Erik, I used refactor fields but with that I'm  not able to get the input value from the user.
I have an input called "Field to Keep" where the user selects the field to be kept. However when I add this input to the Source Expression it understand as a string and all the rows of this column is filled up with the name of the field chosen by the user and not with the actual field value:


texts: QGIS 1.8 or 2.01 select entries in .shp file that match a field entry

texts: I've got a .shp file with 40k points (trees) in it. One field is the species name. Is there a way to select those entries that match a desired field name e.g. 'black walnut' so I can move them to a new layer?


texts: Assigning a polygon shape to density values

texts:
I am pretty new to GIS as a whole. I have a simple flat file in a csv format, as an example:

lat,long,name,val
25,40,abc,10
50,44,def,20

What I am trying to achieve is to assign a 3d polygon in Mapbox such as in this example. While the settings might be quite straight forward in Mapbox where you assign a height and color value based on a data range, it obviously does not work in my case.

I think I am missing out other files such as mentioned in the blog post, like shapefiles or some other file that is required to assign 3d layouts to the 3d extrusion.

I need to know what I am missing out in configuring a 3d polygon, say a cube in Mapbox based on the val data column in my CSV.


texts: Extracting extents of the map view in QGIS

texts: Status bar in QGIS can display min/max of x/y for extents of the map view when panning and zooming. How can these values be extracted and used in expressions? The manual’s expression page doesn’t seem to have variables for current map view extents.


texts: So I have a KML file which has attribute data(Inside the  tag with about 12 data points). This file was generated on ArcGIS. I'm trying to write a PHP script to parse this KML file and insert the datapoints within the ExtendedData tag into a db. How do i do this?


texts: Coordinates of csv file not showing in required location

texts:
I have csv file like this:

I want to plot these places in a map. I have added these coordinates with the help of google earth and tried to plot these coordinates in both QGIS and ArcMap. But, the location is showing in different place than the actual location. How can I solve it?

In QGIS,the x field is longitude and yfield is latitude  and the CRS is EPSG 4326.


texts:
With the following code, I did not succeed to get the first group name of the tree.

from qgis.core import *
from qgis.gui import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *
import os
#coding=system

def group_name():
  root = QgsProject.instance().layerTreeRoot()
  print root.name()

group_name()


texts: How to auto populate multiple fields in attribute table?

texts:
i have Feature Class Buildings, it has many fields like this(Region_id,Region_Name,Governorate_Name,Village_id,Village_Name,Building_id,Building_Hight,Building_Area,Floor_No,Residential_Units_No_Commercial_Units_No,Administrative_Units_No_Other_Units_No)
my challenge is to reduce the time for populate values, i need when i finished create building polygon, then select it to auto populate parameters of attribute values.


texts: GeoServer/OpenLayers: Overlay WMS on Bing Map

texts:
I am coming from the ESRI side and learning GeoServer/OpenLayers. It seems the only way to overlay one layer on Bing or the like is to use WFS. In my opinion, WFS transfers a lot of data to the client side and the client has to render it. It's not that efficient when all I need is to have one layer on the top of the base map. allOverlays=true is too extreme for my case.

Is there a way to overlay it as WMS? Similar to that in ArcGIS Server.

Thanks

Here is my code to try to overlay wms on Bing. The wms layer is not displaying. What's wrong?

        var map;

    function init(){
        var restrictedExtent = new OpenLayers.Bounds(-180, -90,
            180, 90);

        map = new OpenLayers.Map("map");
        //map.allOverlays = true;

        map.addControl(new OpenLayers.Control.LayerSwitcher());

        var shaded = new OpenLayers.Layer.VirtualEarth("Shaded", {
            type: VEMapStyle.Shaded
        });
        var hybrid = new OpenLayers.Layer.VirtualEarth("Hybrid", {
            type: VEMapStyle.Hybrid
        });
        var aerial = new OpenLayers.Layer.VirtualEarth("Aerial", {
            type: VEMapStyle.Aerial
        });

        var states = new OpenLayers.Layer.WMS(
            "States", "http://localhost:8080/geoserver/wms",
            {
                layers: 'state',
                styles: '',
                format: 'image/png',
                transparent: true
            }, {
                opacity: 0.5,
                singleTile: false,
                ratio: 1,
                isBaseLayer: true
            }
        );

        map.addLayers([aerial, hybrid, shaded, states]);

        map.setCenter(new OpenLayers.LonLat(-110, 45), 3);
    }


texts: Has MaxMind's GeoLiteCity IP geocoding database been reduced in accuracy?

texts:
I took a copy of MaxMind's free city-level geocoding database a year or so ago, and the file size was roughly 30MB.

Checking their site again today (July 2012), the equivalent file is now only 17MB, approximately half the size.

I cannot find any reference to this online. Does anyone know if they've reduced the quality or range of their free offering, or whether it's just an improvement in the binary encoding/layout of the data that resulted in the smaller file.


texts: Google Earth Engine: Extract values from FeatureCollection as list

texts:
I'm trying to extract values from a FeatureCollection as a list. Here's some example code that I'm using in Earth Engine's code editor:

var viirs =
    ee.Image('NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG/20120401').select('avg_rad');

var maineCounties =
    ee.FeatureCollection('ft:1S4EB6319wWW2sWQDPhDvmSBIVrD3iEmCLYB7nMM').
    filter(ee.Filter.eq('StateName', 'Maine'));

var maineMeansFeatures = viirs.reduceRegions({
  collection: maineCounties,
  reducer: ee.Reducer.mean(),
  scale: 30
});

The "maineMeansFeatures" FeatureCollection has a column called "mean" -- I trying to extract the values in "mean" as a list (so, here, it'd be a list of 16 values -- one for each feature).

Any ideas on how to do this? (Note: I'm very new to javascript so I'm not that familiar with FeatureCollections).


texts: Embedded quotes in arcpy.CalculateField_management expression

texts:
I'm using arcpy.CalculateField_management to populate a field.  The expression:

"'https://anystate.gov/someportal/DisplayQueryResults.aspx?Prompts=Spill%20Number="!ID!"&Database=OfficeDocs&Query=DEP_Spill_Projects&QuerySet=Portal_Queries&User=Portal.DEP&Password=DEPPortal1'"

will work when entered into the Calculate Field Tool.  The url is completed by pulling in values from the ID field.

However, the Python code snippet pasted into my script results in a syntax error. It's worth noting the expression is bounded by """"'....'"""" (4 double quotes outside a single quote.

I tried Midavalo's suggestion a couple of different ways.
First attempt set

expression = 'https://anystate.gov/someserver/…{0}&Database=OfficeDocs&Query=DEP_Spill_Projects&QuerySet=Portal_Queries&User=Portal.DEP&Password=DEPPortal1'.format(!ID!)

Resulting in invalid syntax
Based on the Python link, the second try included single quotes:

'https://....'.format('!ID!')

But the result was a syntax error


texts: Why is there a difference between NDVI using Sentinel 2, in ArcMap & QGIS?

texts:
I have been trying to calculate the NDVI for a scene using Sentinel 2 bands 4(red) & 8(NIR).
I first did it in ArcMap and then QGIS.
In both softwares, I used the raster calculator and included the basic NDVI formula = (NIR-Red)/(NIR+Red).

ArcMap gave me a range of 0-9221 (this varies with the stretch applied) while QGIS gave me more plausible values from 0.0258-0.5255.

I have read briefly that Sentinel 2 is TOA reflectance and has a scaling factor of 10000, which apparently can be found in the xml file that gets downloaded with the images. I never found this scaling factor mentioned in the xml file I downloaded.

Why is there this difference? How is the calculation being done on QGIS?


texts: UTFGrid grid map generated by TileMill v0.10.1 looks off from the actual map

texts:
I have been trying to create UTFGrid from shapefiles using TileMill v0.10.1. For some reason the grid map it generated looks like it is stretched from the original map.
For Example:

This is the Original Tile:

What TileMill generated:

The top half of the image is correct, however, the bottom half of the character map becomes all '?' and the layout still follows the original image but stretch to the right.

This only happens when the UTFGrid is exported in mbtiles, when using UTFGrid interactivity inside TileMill's GUI, it works perfectly...


texts: Is there a way to traverse and select all vectors connected to a starting vector?

texts:
Preferably using QGis I'd like to remove errors from a road network .shp file and was wondering if there was a away to perform a flood fill or traversal on a vector layer to select all connected vertexes from the starting vertex to do basic validation and remove unconnected roads.

Thanks.


texts: Getting OSM maps in Equirectangular Projection via osm2pgsql

texts:
I'm trying to import some OSM data into psql as equirectangular projection using osm2pgsql with the following command:

osm2pgsql -c -d osm -S /usr/local/share/osm2pgsql/default.style -l /OPENSTREETMAP/south-america-latest.osm.pbf

However, it fails with this error:

Projection code failed to initialise

I have installed the latest version of the PROJ framework from here.

What am I doing wrong?

EDIT 1:

Running 'select * from spatial_ref_sys where srid=4326;' on the database produces this:

 srid | auth_name | auth_srid |                                                                                srtext                                                                                |              proj4text
------+-----------+-----------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+--------------------------------------
 4326 | EPSG      |      4326 | GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]] | +proj=longlat +datum=WGS84 +no_defs
(1 row)


texts: I have this line. You can't see the rest of the line in the picture but It's only one line and I want to connect the two ends to close it and convert it into a polygon. I don't know how to connect it. I'm using QGIS 3.10


texts: Mapbox Studio for macOS renders GeoTiff into webp with black background - suggestions?

texts:
We are converting a 11 MB GeoTiff, 324511630_Descanso_FSTopo.tif,
to post to Mapbox, using Mapbox Studio Classic macOS.

The issue is when the Mapbox tileset is posted in Mapbox.com/studio, there is a black border that extends beyond the original GeoTiffs bounding box.

Why does a black border appear, and why does it extend past the bbox?

The GeoTiff is on the right hand side of this image.  The gray back ground is the Mapbox style light.

Steps to reproduce

wget http://data.fs.usda.gov/geodata/rastergateway/data/32116/fstopo/324511630_Descanso_FSTopo.tif
Mapbox Studio Classic for macOS, New Source > Blank Source > New Layer > Select 324511630_Descanso_FSTopo.tif
Upload to Mapbox Settings > Upload
Alternativly, set the zoom to 11-12 for a faster testing process
The upload is successful with this metadata as reported by Mapbox.com/studio

Format | Type   | Size     | Zoom Extent | Bounds
webp   | raster | 19.28 MB | z11-16      | -180.0,-85.1,180.0,85.1

Notes and Observations

The GeoTiff is georeferenced and is 7.5 minutes, so the bounding box does not match the GeoTiff bounding box
The resulting Mapbox tile is webp, which is what we want as we like the size and performance vs png.

What we've tried

As a comparison, We have modified gdal2tiles.py to use the WEBP tile driver, and it also has the same black border affect.
When rendering with gdal2tiles.py with PNG tile driver, there is no issue with the black ground, but the file size in Mapbox.com/studio is > 50 MB.
from gdalinfo we see this data ColorInterp=Palette & Color Table (RGB with 256 entries) 0: 55,53,53,255
We are aware of the gdal tools, but we don't know how to set the NODATA given the gdalinfo.
gdal_translate -a_nodata 0,0,0,255 in.tif out.tif  # QGIS only accepts one parameters, so this must be wrong

and we tried this, but produces the same issue

gdal_translate -of vrt -expand rgba 324511630_Descanso_FSTopo.tif rgba.vrt
gdal2tilesp.py -z 11 -f WEBP --srcnodata 0,0,0,255 rgba.vrt
mb-util rgba rgba.mbtiles
Upload rgba.maptiles
We also tried -expand rgb instead of rgba.  No difference

gdalinfo 324511630_Descanso_FSTopo.tif

Driver: GTiff/GeoTIFF
Files: 324511630_Descanso_FSTopo.tif
       324511630_Descanso_FSTopo.tif.aux.xml
Size is 9598, 11366
Coordinate System is:
GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS 1980",6378137,298.2572221010002,
            AUTHORITY["EPSG","7019"]],
        AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4269"]]
Origin = (-116.625000000000000,32.875000000000000)
Pixel Size = (0.000013023968607,-0.000010998187728)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=PACKBITS
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-116.6250000,  32.8750000) (116d37'30.00"W, 32d52'30.00"N)
Lower Left  (-116.6250000,  32.7499946) (116d37'30.00"W, 32d44'59.98"N)
Upper Right (-116.4999959,  32.8750000) (116d29'59.99"W, 32d52'30.00"N)
Lower Right (-116.4999959,  32.7499946) (116d29'59.99"W, 32d44'59.98"N)
Center      (-116.5624980,  32.8124973) (116d33'44.99"W, 32d48'44.99"N)
Band 1 Block=9598x1 Type=Byte, ColorInterp=Palette
  Min=0.000 Max=14.000
  Minimum=0.000, Maximum=14.000, Mean=10.365, StdDev=3.715
  Metadata:
    STATISTICS_MAXIMUM=14
    STATISTICS_MEAN=10.365010457425
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=3.7148361607928
  Color Table (RGB with 256 entries)
    0: 55,53,53,255

    1: 238,51,56,255
    2: 253,179,56,255
    3: 187,139,90,255
    4: 132,119,116,255
    5: 207,182,159,255
    6: 187,195,168,255
    7: 253,202,182,255
    8: 234,214,190,255
    9: 200,230,209,255
   10: 253,222,210,255
   11: 220,221,223,255
   12: 0,176,239,255
   13: 201,217,240,255
   14: 255,255,255,255
   15: 0,0,0,255
   16: 0,0,0,255
   17: 0,0,0,255
   <snip>  entries 18 through 254 are identical
   255: 0,0,0,255


texts: How to combine color ramps for two different layers?

texts: I have two layers that contains 5 polygons each with their own values. I would like to make a global color ramp that encompasses both layers. But I am not sure how I can do this. I have tried setting the number of classes under gradient colors to 10 for each layer and create my own color ramp, but ArcGIS automatically resets this back to 5 as there are only 5 values for each layer. The images below display my two layers: 0.25 degrees and 0.5 degrees.


texts: How to programmatically identify and rename duplicate records?

texts:
I'm trying to identify duplicate records from a field and then change these records so they become unique (eg A,B,C, or 1,2,3) and assign this to the end of the text (eg F3241 A56_A or F3241 A56_1).
I've tried the scripts that are in this forum but can't get these to work or not sure if they will give me the result I'm after.
Any input into to this would be greatly appreciated.


texts: When Zooming Personal Geodatabase Feature Class one or more layers failed to draw?

texts: When i joined two different feature class attribute table while i have to zooming the PGDB feature class , drawing error will come. I have attached an example image to show the problem.


texts: OpenLayers function(mergeNewParams) is overwriting previous function call

texts:
I am trying to parse both the values to the corresponding words to fetch data from GeoServer and return me the value, like this:

route.mergeNewParams({viewparams: "'word1:+source;word2:+target"});

But its not working. In fact no values(source and target) are being assigned to the respective words(word1, word2). I have also tried it to do this way:

route.mergeNewParams({viewparams: "word1:"+source});
route.mergeNewParams({viewparams: "word2:"+target});

But in this case last 2nd call of function is overwriting the 1st call.


texts: Opening Esri fGDB Feature Class in MapInfo Professional 11.0.3

texts:
I would like to open ArcGIS file Geodatabase feature class and other such fGDB objects with MapInfo Professional. NB: converting to Shapefiles or personal mdb is not an acceptable solution. I need it to be seamless and easy for users, just like the ArcSDE and .mdb methods (pictured below).

Any ideas are welcome


texts: Georeferencing raster using GDAL and Python?

texts:
I want to georeference a raster using python and GDAL. My current approach is to call gdal_translate and gdalwarp using os.system and an ugly list of ground control points. I'd really like a way to do this natively within python.

This is the current process I am using:

import os

os.system('gdal_translate -of GTiff -gcp 1251.92 414.538 -7.9164e+06 5.21094e+06 -gcp 865.827 107.699 -7.91651e+06 5.21104e+06  "inraster.tif" "outraster1.tif"')
os.system('gdalwarp -r bilinear -tps -co COMPRESS=NONE "outraster2.tif" "outraster3.tif"')

There is a previous question and answer from 2012 which states that gdal_translate can be accessed after importing gdal. I'm not sure if is obsolete, or whether it is wrong but when I run from osgeo import gdal I do not see gdal.gdal_translate as an option.

I don't know if it exists but I would love if I could translate and reproject rasters in a pythonic way. For example:

# translate
gcp_points = [(1251.92, 414.538), (-7.9164e+06, 5.21094e+06)]
gdal.gdal_translate(in_raster, gcp_points, out_raster1)

# warp
gdal.gdalwarp(out_raster1, out_raster2, 'bilinear', args*)

Is such an approach possible?


texts: Most useful program to create an accurate city map for printing?

texts:
I googled around and several sources said this would be a good place to ask about mapping so here goes. Thanks in advance for reading!

I want to make my own city map (specifically Paris) as a gift to the girlfriend before she visits later this year. I would like for it to be aesthetically pleasing and as professional as possible since I'm planning to make a foldout map-sized professional print of the final product and perhaps place it inside a notebook containing other useful travel things.

At first, I thought of hand-drawing the map, but figured it would be much too time-consuming to draw every street in Paris by hand. Then I considered using Illustrator or a similar image creation/editing program to fill in the details of an available map. But, I think using some sort of more professional mapmaking software with GIS capabilities would actually be best, since I'd like to accurately place landmarks, sights, museums, and other destinations of interest, as well as metro stations, street names, and the other bread and butter features of a good city map.

I know I could easily find a good city map of Paris, but I'd like to make it a more personal gift, as well as do my own customizations, like specifying the color scheme and selecting which features to include and label. Also, I really enjoy geography and love to learn new things, so why not get my hands dirty with GIS and mapping while making this gift?

So, what are your suggestions for a good program to use? One which allows accurate placement of features and landmarks, has color and style palettes and tools to customize the appearance, and exports to high quality image files for large-format printing.

Thanks for your input!


texts: Using Python to convert line shapefile to raster, value=total length of lines within cell

texts:
I have already read the similar question Convert line shapefile to raster, value=total length of lines within cell using R.
And I want to realize the function using Python.

===== My attempt ====

generate the grid network

 ### (1) 140 rows x 187 columns (1km x 1km grid)
 lon_x = np.linspace(xc1,xc2,187)
 lat_y = np.linspace(yc1,yc2,140)

 ### (2) empty array for saving the length data within the grid
 sh = (140*187,2)
 grids = np.zeros(140*187*2).reshape(*sh)
 grids_road = np.zeros(140*187)

 ### (3) generate the networks and k-d tree structure
 xx = lon_x
 yy = lat_y

 k = 0
 for j in range(0,yy.shape[0],1):
     for i in range(0,xx.shape[0],1):
         grids[k] = np.array([xx[i],yy[j]])
         k+=1

####  the k-d tree algorithm I used can be applied for finding the start/end grid agree with one road
T = spatial.KDTree(grids)

Loop all my road polyline => intersection with grid => saving the length

### define the search radius
x_delta = (lon_x[24] - lon_x[23])
y_delta = (lat_y[24] - lat_y[23])
R = np.sqrt(x_delta**2 + y_delta**2)

### loop the road(I save the polyline as pandas dataframe which contain the start point x1,y1 and end point x2, y2 and length)

for i in range(0,len(road),1):

    sta = sorted(T.query_ball_point([road.x1.iloc[i],road.y1.iloc[i]],r=R))[0]
    end = sorted(T.query_ball_point([road.x2.iloc[i],road.y2.iloc[i]],r=R))[0]

    dt = (grids[end][0] - grids[sta][0])/x_delta
    dx = round((grids[end][0] - grids[sta][0])/x_delta)
    dy = round((grids[end][1] - grids[sta][1])/y_delta)

    # using shapely to change the road into line shapefile
    line = [(road.x1.iloc[i], road.y1.iloc[i]), (road.x2.iloc[i],road.y2.iloc[i])]
    shapely_line = shapely.geometry.LineString(line)

    # if the road extend just in one grid, there'll be no intersection
    if road.distance.iloc[i] < 1.0:
         index = sta
         grids_road[index] = road.distance.iloc[i]
    if road.distance.iloc[i] >= 1.0:
       if (dx > 0) & (dy == 0):
          for j in range(0, int(dx),1):
              k = 0
              x1,x2,x3,x4 = grids[sta][0] + j*x_delta, grids[sta][0] + (j+1)*x_delta,\
                       grids[sta][0] + (j+1)*x_delta, grids[sta][0] + j*x_delta
              y1,y2,y3,y4 = grids[sta][1] + k*y_delta, grids[sta][1] + k*y_delta,\
                       grids[sta][1] + (k+1)*y_delta, grids[sta][1] + (k+1)*y_delta,
              polygon = [(x1,y1),(x2,y2), (x3,y3), (x4,y4), (x1,y1)]
              shapely_poly = shapely.geometry.Polygon(polygon)
              intersection = list(shapely_poly.intersection(shapely_line).coords)
              ds1,ds2 = intersection[0],intersection[1]
              index = sta + j
              grids_road[index] = vincenty(ds1,ds2).miles

Then I classify the road into different scene based on the trends(dx, dy can represent this property). But I found this solution is too rigid, and I want to achieve the target with more  efficient  code.

Wish for your guide!

The figure shows my early progress which contain some  vertical and horizontal road.
http://i4.tietuku.com/028bdd8ef1207557.png

Add

This method was too primary comparing with the method using R.
And, when the road wasn't parallel to Lon or Lat, the code using this method was too hard to write for different situations.
So, I was think whether python has some similar functions or not?


texts: Geometric Calculations for Group of Pixels

texts:
I have group of pixels that I created by region group tool and I would like to measure the length and the width of each pixel group. I need max length and the width at half of the length of each group

My main goal is to extract the structure of each pixel group i.e whether it is linear or other pattern.

Would you recommend any tool for that?


texts: Calculating shortest distance between two points on water, avoiding land?

texts:
Can someone explain to me (in beginner's terms) how to calculate the distance between two points that are on water, avoiding land?

For example, these points are inshore in a bay, so I don't want the path to go over any of the bends in the coastline.


texts: OpenLayers not requesting tiles at previously visited zoom levels

texts:
I have a web app with Maps on several pages using OpenLayers as a client and GeoServer as a backend. The layers are WMS tiles.

On one page, if you visit it, and go to several zoom levels (of the WMS Tile layer),

and then say go into GeoServer and change the styling on that layer.

Then come back to this page, the style will not change for these zoom levels, and watching the browser requests, requests do not occurr at previously visited zoom levels.

I have Response Cached Headers OFF, Tile Caching OFF on the Layer.

Why is the request not being sent?  Note the web pages are asp.net, requested via AJAX.


texts: Creating raster stack from xarray in Python

texts:
I am trying to create a raster stack from an xarray dataset which I obtained from multiple NetCDF files.
import netCDF4 as nc
import rasterio as rio
import numpy as np
import xarray as xr
import os

fpath = '/home/sst/2015'
pattern = '*.nc'
filelist = []

for root, dirs, files in os.walk(fpath):
for name in files:
    filelist.append(os.path.join(fpath,name))

ds = xr.open_mfdataset(eflist, concat_dim='time', parallel = True) # netCDF data

ds_data = ds.sel(time='2015')['SST']   # xarray dataset with dimension 365x3600x7200.

By looping over the 3D xarray following code writes 365 files to HDD and later I can read and stack those files.
from rasterio.transform import from_origin
transform = from_origin(-180,90, 0.05, 0.05)

fpath = '/home/sst/sst_tif'
fname = 'sst_array_'
extname = '.tiff'
timedf = ds.time    # time dimension to loop over

for i in range(len(timedf)):
    np_array = np.array(ds_data[i])
    iname = str(i)
    fwname = fpath + fname + iname + extname

    sst_tif = rio.open(fwname,
                    'w',
                   driver = 'GTiff',
                   height = 3600,
                   width = 7200,
                   dtype = np_array.dtype,
                   count = 1,
                   crs = 'EPSG:4326',
                   transform = transform)

    sst_tif.write(np_array, 1)
    sst_tif.close()

However this takes a very long time as it creates 365 files on HDD. Is there any way I can create this raster stack in memory and do further processing without having to write it to file/s on HDD? I am trying to obtain functionality similar to that of stack available in raster package of R.


texts:
I have a basic tracking system i've written with OpenLayers using google and openstreetmap and I want to use geocoding to get street names etc

However, I'm a bit confused.
I was hoping there was a easy/free way with OL and OSM however, i'm not sure there is

I've tried using Google, but it doesn't work with bulk data.
It allows me to get about 3 or 4 back and then doesn't work for the rest..... i think they require a wait period between each request

Can anyone recommend a free way to geocode without Google?
please be patient with my lack on knowledge in this areas


texts: how to change the Geometry of Polyline in QGIS using Python

texts:
I have two layers with multiple features, BlackLayer and YellowLayer. As i mentioned in the above image there are one yellow feature on two black features.I would like to change the geometry of yellow feature to black feature as i shown in below image.

what i tried is:

lineLayer = iface.activeLayer()
provider = lineLayer.dataProvider()

spIndex = QgsSpatialIndex() #create spatial index object

feat = QgsFeature()
fit = provider.getFeatures() #gets all features in layer

#insert features to index
while fit.nextFeature(feat):
    spIndex.insertFeature(feat)

pt = QgsPoint(-0.00062201,0.00001746)

#QgsSpatialIndex.nearestNeighbor (QgsPoint point, int neighbors)
nearestIds = spIndex.nearestNeighbor(pt,1) # we need only one neighbour

featureId = nearestIds[0]
fit2 = lineLayer.getFeatures(QgsFeatureRequest().setFilterFid(featureId))
ftr = QgsFeature()
fit2.nextFeature(ftr)
vertices = ftr.geometry().asPolyline()
polyline = QgsGeometry.fromPolyline(vertices )
provider.changeGeometryValues({ftr.id(): polyline})

Please help me with this ..!!!


texts: Creating subset for doing Co-Kriging interpolation and Crossvalidation?

texts:
I have a data set with only 47 soundings, so a small amount of data available. With this information I am trying to produce a raster surface through interpolation. Unfortunately, a data analysis showed that my data is not autocorrelated. So I chose a multivariate approach with the Geostatistical Analyst: Co-Kriging.

Now I am wondering if it is necessary to create a subset before doing the interpolation, because as far as I understand, the crossvalidation which can be done with the Geostatistical Analyst applies the leave-one-out technique. This means I can use the whole data set and still validate the results appropriately, right?


texts: Snap LineStrings considering direction in QGIS 3

texts:
I have 2 projected shape files: Roads (red) and public transit routes (blue):

I want to snap the PT routes to the roads, but because of inaccurate marks sometimes the routes are snapped to the wrong road segment.
I wish to minimize the errors by considering the direction while snapping, and not only distance. Is there any way to do that in QGIS?


texts: I'd like to take a test about the connection to the database in QField in my smartphone Android but it failed; so I cannot view the layer Postgis on the map. The database Postgresql is installed and connected to my laptop, SSL is disabled, the QGIS project has the DB credentials stored and I exported the project via Qfield sync as a package for QField (setting "no action" for the PostGIS). Any suggestions?


texts: how to change marker angle in openlayers?

texts:
I want show one marker with normal angle 90 degrees(it is default of openlayers) and other marker with 120 degrees in openlayers. To differentiate with the normal ones.
how can i do this.


texts: How to add multiple Landsat image on GEE?

texts:
I want to use Landsat image from GEE with area specified within an input. But when I run the script it only shows a small part of the area (the intended area of interest is the one in blue). The code is as the following
var countries = Indonesia

var Nama_Provinsi = ['Kalimantan Selatan']
var geometry = countries.filter(ee.Filter.inList('NAME_1', Nama_Provinsi));
var area1 = ee.FeatureCollection(geometry);
Map.centerObject(geometry, 8);
Map.addLayer(geometry,{color:"Blue" },"Provinsi");

var L8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var RGB_vis = {min: 0, max: 0.3, bands: ['B4', 'B3', 'B2']};
var filtered = L8
.filterDate('2020-01-01', '2020-12-31')
.filterBounds(geometry)
.sort('CLOUD_COVER', true);
var image = ee.Image(filtered.first());
Map.addLayer(image, RGB_vis,'RGB', false );

// Compute the Enhanced Vegetation Index (EVI).
var nir = image.select('B5');
var red = image.select('B4');
var blue = image.select('B2');
var evi = image.expression(
  '2.5*(B5-B4)/(B5+6*B4-7.5*B2+1)', {
    'B5' : nir,
    'B4' : red,
    'B2' : blue
  });
// Display the result.
Map.centerObject(image, 8);
var eviPara ={min: -1, max: 1, palette: ['blue', 'white', 'green']};
Map.addLayer(evi, eviPara, 'EVI');

And the result I got is as the picture.


texts: Landsat 8 extracted pixel values for up North: 'pixel_qa' value of 322

texts: I am trying to get landsat reflectance values corresponding to field plots spread all over Finland using the Google Earth Engine. As you know, Finland is up North, between latitudes of ~60 and 70 degree N. I am extracting from this landsat collection. Upon examining the extracted values, I noticed that most (99.9%) of the pixels have a 'pixel_qa' value of 322. From my understanding of the bitmask, I gathered that the value of '322' means that the 'cloud confidence' and the 'cirrus confidence' was low. Why is this? Can I still assume that the cloud/noncloud classification is good enough?


texts: How can I access the user's location in my Google Earth Engine script?

texts:
I'm working in Google Earth Engine and I want to access the user's location to then be able to use the position when visualizing data layers. When running the code below, GEE doesn't seem to be able to understand the navigator which contains the geolocation info I'm looking for. Same goes for document. This makes me wonder whether GEE even works in JS since it doesn't understand the queries(?) Anyone with a solution to access the user location?

var locationButton = ui.Button('User location');
locationButton.onClick(function() {

  var startPos;
  var geoSuccess = function(position) {
    startPos = position;
    document.getElementById('startLat').innerHTML = startPos.coords.latitude;
    document.getElementById('startLon').innerHTML = startPos.coords.longitude;
     if (navigator.geolocation) {
      console.log('Geolocation is supported!');
    }
    else {
      console.log('Geolocation is not supported for this Browser/OS.');
    }
  };
  navigator.geolocation.getCurrentPosition(geoSuccess);

});
Map.add(locationButton);


texts:
QGIS does not recognize .prj file, which is part of a shapefile.

ArcMap does automatically recognize .prj file, when you add shapefile to ArcMap layer, you click property, coordinate tab, you can see the correct CRS (coordinate reference system) already set there.

However, QGIS failed to recognize .prj file, instead you have to manually set CRS in QGIS.

But how?


texts:
I am relatively new to GIS. I am trying to compare a DEM that I developed through photos from an aerial drone and a DEM that we developed through traditional surveying methods.

Exactly how would I get the RMSE? I saw some similarly asked questions that explained how to use R to get it, but I am very unfamiliar with any sort of coding. Is there anyway to get the RMSE through Spatial Analyst or other already available tools?


texts: How to make zooming for OpenStreetMap in OpenLayers more sensitive to mouse wheel scrolling?

texts:
I have been testing OpenLayers for a while and I notice that the scrolling on OpenStreetMap layer does not make the map zoom in or out as much as I would like it to be. On most occasions, I am only able to zoom in or out by only one level at a time, regardless of the scroll clicks. However, when the base layer is switched to Google Maps, I am able to zoom in and out more with more scroll clicks. You can see the difference with this example.

At first, I thought this could be due to some inherent problem with OpenStreetMaps, but after testing its own map on the main site, there seem to be no problem with zooming through multiple levels with a mouse scroll. Am I missing something in the OpenLayers control setting to enable such responsive scrolling behaviour?


texts: Python and GDAL: NoneType error while reading shapefile

texts:
I am a newbie to GDAL. I have a shapefile of Britain which I want to rasterize. This file can be obtained from here selecting United Kingdom as Country and Administrative areas as Subject.

Following the GDAL tutorial ("Convert vector layer to array"), I have this script:

import pyproj
import osgeo
import numpy as np
from osgeo import gdal
import os
import ogr, gdal
from gdalconst import *

#The dir and file to work with
mydir=r'C:\Users\GBR_adm'
os.chdir(mydir)
vector_fn = 'GBR_adm0.shp'

# Define pixel_size and NoData value of new raster
pixel_size = 25
NoData_value = 255

# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
source_srs = source_layer.GetSpatialRef()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('MEM').Create('', x_res, y_res, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])

# Read as array
array = band.ReadAsArray()
print array

This script raises an error: AttributeError: 'NoneType' object has no attribute 'SetGeoTransform', and this relates to type(target_ds) being NoneType. I have read somewhere that this means that GDAL has failed to load the shapefile. Why does this happen, and how to avoid it?


texts: Visualizing two different datasets amicably when they are measuring two different things

texts:
I am working on a story map in AGOL that is comparing the percentage of a population that lives in a house that has been deemed transit dependent versus the concentration of a population that lives in a house that has been deemed transit dependent. Each of these are done by census tract.

The purpose of the story map is to show that while there are areas that have been deemed by the area density map as not as significant because of the overall size of the census tract there are areas that have significant portions that are not occupied by residents, but the overall percentage is high.

How do I ensure that these two scales are evenly distributed while at the same time retaining the simplicity of user interpretation (not using a standard deviation scale)

The below image is the scale that I am currently working with that measures percentage of the population versus residents by square mile.


texts: Exporting vector and raster in native projection

texts:
I'm trying to export raster and vector in native projection and resolution. And I have made a code for it. After exporting both raster and vector versions of data, I have added them to ArcMap, but raster version does not lie geographically correct, it is shifted. I have tried to find a solution and I could not. Anyone you know why it is happening and how to solve it?

Please find the link below:
https://code.earthengine.google.com/33a783d191fefea8ef6c5956ddb01a01

var modis_col = ee.ImageCollection('MODIS/006/MOD13Q1').select('NDVI')
    .filterDate('2000-04-01', '2000-06-10');

var modis_ndvi = modis_col.max().divide(10000).gte(0.4).reproject(modis_col.first().projection());

var vis = {min: 0, max: 10000, palette: [
  'FFFFFF', 'CE7E45', 'FCD163', '66A000', '207401',
  '056201', '004C00', '023B01', '012E01', '011301'
]};

Map.addLayer(modis_col.first().clip(roi), vis, 'NDVI');
Map.addLayer(modis_ndvi.clip(roi));

var modis_ndviMasked = modis_ndvi.updateMask(modis_ndvi.gte(0.5));

var modis_ndvi_vector = modis_ndviMasked.reduceToVectors({
  geometry: roi,
  crs: modis_col.first().projection(),
  geometryType: 'polygon'
});

Map.addLayer(modis_ndvi_vector)

Export.table.toDrive({
  collection: modis_ndvi_vector,
  description: 'modis_ndvi_classed_2000_shp',
  fileFormat: 'SHP',
});

Export.image.toDrive({
  image: modis_ndviMasked,
  description: 'modis_ndvi_classed_2000_tif',
  crs: modis_col.first().projection().getInfo().crs,
  crsTransform: modis_col.first().projection().getInfo().transform,
  region: roi
});


texts: Select points approx. 2000 metres from another point along a river?

texts:
I have a series of points representing a river network (See Image Above). Each point contains attributes on Latitude, Longitude, Z value and Length calculated using flow accumulation. The river network was created from Digital Elevation Models using hydrology tools in ArcGIS (Fill, Flow Direction, Flow Accumulation, Strahler Order > 4). I also have a polyline shapefile that also runs along this river network.

Starting at the first point at the beginning of each network I would like to determine another point "downstream" that is approximately 2000 metres away (Or closest to 2000 metres). I would then like to calculate channel slope based on the difference in elevation divided by the difference in length (2000m) using the two points.

I am aware that at the junctions there will be an issue, any suggestions on accomplishing this task in GIS to determine the point 2000m downstream and continuing downstream at junctions?

I am using ArcGIS 10.1 with an ArcInfo license and also have access to QGIS and GrassGIS. I am familiar with using python scripts as well.

Edit: Progress

The stream polyline has the following attributes:

x1: x coordinate of start point (Double)

y1: y coordinate of start point (Double)

x2: x coordinate of end point (Double)

y2: y coordinate of end point (Double)

length: length of polyline

shape: shape of polyline

startxy: [x1]&", "&[y1] (Text)

endxy: [x2]&", "&[y2] (Text)

I am using the tool to make arrays

start = arcpy.da.FeatureClassToNumPyArray("polyline",["x1","y1"])

Dictionary

G ={}
>>> for (st,sh,le,en) in zip (start,shape,length,end):
...     G[st]=(st,le,en)

Algorithm

def Downstream (p, x, G):
...     e = G[p]
...     if (IsNull(e)):
...         return ("Not Found")
...         if (x < 0):
...             return ("Invalid")
...             if (x < e.length):
...                 return (Along (e.shape, x))
...                 return (Downstream (e.end, x-e.length,G))

The Downstream function is receiving a runtime error where there is a key error.
I have used G.keys() to see the keys, but they are not working when inputted in the function.


texts: Multimodal network v.s. single mode network

texts:
I created a multimodal network with rail, bus, and pedestrian network, and generated a 15 min service area (blue in the attached figure). I also created a single mode network with only pedestrian network, and got another 15 min service area (yellow). When I compared them, I found some area that used to be covered by 15-min walk has not been covered by the multimodal system, such as the north branch of the yellow. The issue gets really problematic as the gap increases when 30 mins and 45 mins are calculated.

Another problem is that the shape of multimodal 15-min service area is very segmented and does not make sense.

My network is clean that all road/bus route/rail route are made of links and I selected honor as the connectivity policy. No traffic restrictions applied yet at this moment.


texts: Calling gdal_pansharpen.py from another Python script?

texts:
Has anyone been able to successfully call gdal_pansharpen.py from your own Python script?

I have successfully used rasterio to combine my R, G, and B bands, but would like to perform a pansharpen on the final image using the gdal_pansharpen command.

How did you manage to make the call?


texts:
I have an SRTM-derived hillshade (30x30m) and a Quickbird scene (0.6x0.6m).
Is it possible to kind of pan-sharp the hillshade to have the detail but also have the topography?


texts: QGIS 3.12 raster reprojection giving no result

texts:
I would like to reproject my raster DEM from the EPSG:4376 (WGS-84) to the EPSG:27700 (British National Grid 1936). When I use the option: Raster -> Translate (COnvert format), as shown below:

Everything proceeds alright, which finishes the comment:
"Algorithm 'Translate (convert format)' finished"

Unfortunately, it's not translated into the visual result, because the layer is completely transparent.

Further saving the layer as a .TIFF format doesn't help.

What am I doing wrong here then?


texts: Inputting integer values into SQL statement using ArcPy?

texts:
I am defining a function for use in a ArcGIS tool that will verify attributes, catch errors, and obtain user input to rectify those error. I want the tool to select and zoom to the segment that is being currently assessed so that they can make an informed decision. This is what I have been using, and it works well. But the CONVWGID is the variable that will be changing, and I'm not sure how to input that variable into an SQL statement without causing errors.

This is how I had tested the logic:

def selectzoom():
    arcpy.SelectLayerByAttribute_management(Convwks, "NEW_SELECTION", " [CONVWGID] = 10000001")

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

Then I needed to work the variable into the function in order to accept different CONVWGID values, which gives me a Runtime/TypeError that I should have known would happen.

Runtime error -
Traceback (most recent call last): - File "string", line 1, in module - TypeError: cannot concatenate 'str' and 'int' objects

def selectzoom(convwkgid):
    delimfield = '" [CONVWGID] = ' + convwkgid + ' "'
    arcpy.SelectLayerByAttribute_management(Convwks, "NEW_SELECTION", delimfield)

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

And when I alter the delimfield line to change the integer into a string, it selects all of the attributes in the entire feature class. Not just the one that had been passed via the function call.

delimfield = '"[CONVWGID] = ' + str(convwkgid) + '"'

I'm not amazing with SQL and maybe I'm missing something basic with this statement, but I can't figure out why it won't work when I'm basically giving it the same information:

"[CONVWGID] = 10000001"
'"[CONVWGID] = ' + str(convwkgid) + '"'


texts:
i'm struggling to run batch processing multilevel b-spline interpolation for my surface area the below error message keeps on poping up after the process>

Processing algorithm 1/3...
Algorithm Multilevel b-spline interpolation starting...
grid_spline "Multilevel B-Spline Interpolation" -TARGET_DEFINITION 0 -SHAPES "C:/Users/mlamli.v/Downloads/Trial3_B1_1_Site.csv.shp" -FIELD "2011-07-08 " -METHOD 0 -EPSILON 0.0001 -TARGET_USER_XMIN 19.2206320676 -TARGET_USER_XMAX 19.2230436656 -TARGET_USER_YMIN -34.0070661627 -TARGET_USER_YMAX -34.0056905999 -TARGET_USER_SIZE 1e-06 -TARGET_USER_FITS 0 -TARGET_OUT_GRID "C:\Users\mlamli.v\AppData\Local\Temp\processing2e2280fd77be44ffa848faacde5e36a1\78ac882c583747e4a712f0ca8379c16d\TARGETOUTGRID.sdat"

C:\PROGRA~1\QGIS2~1.18\bin>set SAGA=C:/PROGRA~1/QGIS2~1.18/apps\saga-ltr

C:\PROGRA~1\QGIS2~1.18\bin>set SAGA_MLB=C:/PROGRA~1/QGIS2~1.18/apps\saga-ltr\modules

C:\PROGRA~1\QGIS2~1.18\bin>PATH=C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg\shapely\DLLs;C:\PROGRA~1\QGIS2~1.18\apps\Python27\DLLs;C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\numpy\core;C:\PROGRA~1\QGIS2~1.18\apps\qgis\bin;C:\PROGRA~1\QGIS2~1.18\apps\Python27\Scripts;C:\PROGRA~1\QGIS2~1.18\bin;C:\Windows\system32;C:\Windows;C:\Windows\system32\WBem;;C:/PROGRA~1/QGIS2~1.18/apps\saga-ltr;C:/PROGRA~1/QGIS2~1.18/apps\saga-ltr\modules

C:\PROGRA~1\QGIS2~1.18\bin>saga_cmd grid_spline "Multilevel B-Spline Interpolation" -TARGET_DEFINITION 0 -SHAPES "C:/Users/mlamli.v/Downloads/Trial3_B1_1_Site.csv.shp" -FIELD "2011-07-08
____________________________

##### ## ##### ##
### ### ## ###
### # ## ## #### # ##
### ##### ## # #####
##### # ## ##### # ##
____________________________

SAGA Version: 2.3.2 (64 bit)

____________________________
library path: C:\PROGRA~1\QGIS2~1.18\apps\saga-ltr\modules\
library name: grid_spline
library : Spline Interpolation
tool : Multilevel B-Spline Interpolation
author : O.Conrad (c) 2006
processors : 4 [4]
____________________________

Load shapes: C:/Users/mlamli.v/Downloads/Trial3_B1_1_Site.csv.shp...

Parameters

Points: Trial3_B1_1_Site.csv
Attribute: 2011-07-08
Target Grid System: user defined
Cellsize: 0.000028
Left: 19.221356
Right: 19.222112
Bottom: -34.006392
Top: -34.005944
Fit: nodes
Method: with B-spline refinement
Threshold Error: 0.000100
Maximum Level: 11

level:1, error:17, max:46.947204, mean:24.872109

level:2, error:17, max:41.319060, mean:21.044483

level:3, error:17, max:31.586489, mean:13.923085

level:4, error:17, max:22.824775, mean:6.834893

level:5, error:17, max:9.743622, mean:2.226048

level:6, error:8, max:1.612531, mean:0.710186

level:7, error:2, max:0.002017, mean:0.001562

level:8, error:0, max:0.000000, mean:0.000000

C:\PROGRA~1\QGIS2~1.18\bin>" -METHOD 0 -EPSILON 0.0001 -TARGET_USER_XMIN 19.2206320676 -TARGET_USER_XMAX 19.2230436656 -TARGET_USER_YMIN -34.0070661627 -TARGET_USER_YMAX -34.0056905999 -TARGET_USER_SIZE 1e-06 -TARGET_USER_FITS 0 -TARGET_OUT_GRID "C:\Users\mlamli.v\AppData\Local\Temp\processing2e2280fd77be44ffa848faacde5e36a1\78ac882c583747e4a712f0ca8379c16d\TARGETOUTGRID.sdat"
The input line is too long.

C:\PROGRA~1\QGIS2~1.18\bin>exit
Converting outputs
'NoneType' object has no attribute 'crs' See log for more details

Can anyone assist me?


texts: Creating ant task to build ArcGIS Viewer for Flex?

texts:
Has anybody created ant task to build ArcGIS for Flex Viewer?

I want to put ArcGIS Flex Viewer into our Continuous Integration system (Hudson).

So far it's a very difficult task to do.


texts: Count number of pixel for each polygon and have the data in one featurecollection

texts:
I have raster with NDVI values. I have calculate for each plot the number of pixels in different values (no. of poixels lower than 0.2, higher than 0.2 and lower than 0.4 ect....)  based on masks I have created.
My problem is that now I compute it for the polygons and create each time new featurecolelction that contains this count data.
I would like to have one featurecollection that has the data of the counting for each mask layer.

Here you can see my code and try it, the raster and the polygon suppose to be available.

https://code.earthengine.google.com/675b9c9329db3c8e1d556ee577443173

My end goal: To have one Featurecollection which will contain serveral columns, one column on each "count" process with the " maineMeansFeatures" calclations.


texts:
I work with the latest QGIS program on bike trail maps. I monitor my CPU and memory usage and for years the CPU usage hovers around 20%, never to peak or even hit 50%.
This just seems odd that I have to wait for tasks to get done when there is still so much headroom to get it done faster. Is this normal?
I work on a Win 10 computer with an i7, 16MB Ram, new video card etc. it flies.
And yes I have set it to use all cores, work with Geopackage files, clip my maps smaller...
Would setting up a DB with QGIS tap into more power?
The time lost in slow processing and rendering doing 66 maps is a real problem.


texts: Is it possible to lock the layers using Qgis API Python?

texts: I would need to load and "lock" two layers on the map. Locked means that layers cannot be removed from the layer registry. When I say removed I mean removed by the user using the QGIS interface (the remove button). And possibly make them as read only? Only if this would not interfere with the plug-in I am making, which manages and edit layer geometries, through customized GUIs I made.


texts: Returning shape.area in code block of ArcGIS Field Calculator with if statement?

texts:
Using ArcGIS 10.1 for Desktop, I try to calculate a field using an 'if' statement in the field calculator.
Here is my code  :

"
#Expression
MyField_2 = def Reclass(!MyField_1!)

#Parser
Python 9.3

#Code Block
def Reclass(MyField1) :
    if MyField1 == -1 :
        return 0
    else :
        return !shape.area!
"

But I have an error returning !shape.area!.
How can I fix that ?


texts: Calculate NDVI image from R and NIR channels and rescale it to 0-255

texts:
I have images of plants taken in Red(R) and Near End Infrared(NIR) field,

I calculate the NDVI image per pixel using the below formula,

      NDVI = (R-NIR)/(R+NIR) , for all pixels,

I get a range from [-1,1].

Here NDVI image is used to separate soil from plant.

Now How do I convert it to a range from 0 to 255(unsigned int) for image processing.

Any help ?

EDIT

I am coding in C++ and I want a formula that does the conversion so that I can code it.I use OpenCV for image processing.
I am not using any tool to the conversion.


texts:
I am displaying multiple feature types using one layer in OpenLayers, like here:

 var **vector** = new OpenLayers.Layer.Vector("My Layer", {
  maxExtent: new OpenLayers.Bounds(-1.09863,44.42907, 91.01074, 2.05840)
  sphericalMercator: true,
  protocol:  new OpenLayers.Protocol.WFS({
  url: 'http://mydomain/geoserver/wfs',
  featureType: **["mylayer_1","mylayer_2","mylayer_3"],**
  featureNS: "http://mapmap.org",
  srsName: "EPSG:3857",
  version: "1.1.0",
  extractAttributes: true,
  isBaseLayer: false,
  visibility: true
   })
  });

now I want to make layer_2 invisible by setVisibility(false).. but I don't know how to select it. it seems like I just can control the whole layer (vector).

is there any way to control visibility of these layers ?


texts: QGIS: habitat suitability/combining raster layers

texts:
I have created a number of raster layers for different environmental factors (slope, roads, rivers, forest cover) and want to combine them to find areas where the 'ideal' aspects all overlap. For example, where the slope is not too steep, rivers are within 3km, etc.

I have already added buffers to the roads and streams, as well as reclassified the slope and vegetation layers to have categories. All layers are rasters.


texts:
I am creating a map for disaster management -currently Puerto Rico and we need to make it accessible to a team of 100's of people but not make it public. I know we can use the Organisation Tab and add users but this only allows about 20.

Can we just have the web app to ask for a generic login/password that we share with those that need access or ideally link to a social media account that we assign access to?

I am new to ArcGIS Online.


texts:
WMS can be added to Folium Maps. Similarly, is there a way to add ArcGIS MapServer in Folium Map as a layer? For example, I want to add this server to the Folium map
https://ags02.sec.usace.army.mil/server/rest/services/NLD2_PUBLIC/MapServer/7
I looked at the query that QGIS is sending to get the image and it is as follows:

https://ags02.sec.usace.army.mil/server/rest/services/NLD2_PUBLIC/MapServer/export?bbox=-97.163384,28.215345,-96.425103,28.942482&size=795,783&format=PNG32&layers=show:7&transparent=true&f=image

It is not a tile server, so it cannot be added as tiles.


texts: Why is incorrect Travel Time: < 1 minute being reported by ArcGIS Network Analyst Find Route?

texts:
I have configured a network dataset using ArcGIS 10.4.1 for Desktop and am testing it using the Find Route tool.  Results are as expected except that on the Directions tab the Travel Time is always being reported as Travel Time: < 1 minute.

The example below travels along a single poly line which Identify reports has LENGTH = 3061.1 (metres), SPEED = 50 (kmh) and Minutes = 3.67479.  Consequently, I would expect that Travel Time reported would be 3.67479 minutes.

The network dataset properties which I think are relevant are as below and the network dataset has been re-built and ArcMap re-started since I calculated the Minutes field.

Is there an obvious setting/step that I am likely to be missing?


texts: Export the contour lines of the relief from OSM

texts:
My goal is to display the contour lines of the relief in Google Earth.
So, I would like to export from OpenStreetMap the contour lines of the relief of a region.
My goal is to download the data as a kml, to load them into Google Map.
I was hoping a query like that to work in OverPass:
node
  ["contour"="elevation"]
  (42,1,43,1.5);
out body;

It does not. I could not find a more appropriate tag. According to OSM wiki, it does not exist.
Is there any other way?


texts: Create KML from CSV file, with 2 sets of lon-lats (instead of 1)

texts:
Having the dataframe:

lon-1 lat-1 lon-2 lat-2 some_val

And using the R code:

coordinates(mydata)<-~lon_start+lat_start+lon_end+lat_end
proj4string(mydata) <- CRS("+init=epsg:28992")
mydata_ll <- spTransform(mydata, CRS("+proj=longlat +datum=WGS84"))
writeOGR(wwverkeer_ll["some_val"], "mydata.kml", layer="some_val", driver="KML")

I get an unexpected output (kml file) containing:

<Placemark>
  <ExtendedData><SchemaData schemaUrl="#some_val">
      <SimpleData name="some_val">3228</SimpleData>
  </SchemaData></ExtendedData>
   <Point>
<coordinates>3.xxxxxx,47.yyyyyyy,54.zzzzzzzz</coordinates></Point>
   </Placemark>
   <Placemark>

Where I would actually like to have 4 coordinate points, instead of 3.
And some certainty on this being the correct coordinates.


texts: Do QGIS Field Calculator xat(n) and yat(n) functions work in version 1.7.x?

texts:
I'm trying to use the xat(n) and yat(n) in the Field Calculator but can't get them to work.
Does anybody have an example how to use them?
I need the start and end coordinates for line segments, does anybody know how to get those?

In the below list the xat() and yat() functions are listed.
http://hub.qgis.org/wiki/17/List_of_Field_Calculator_Functions


texts: Qt/C++ Hello World application with QGIS not compiled: qgsconfig.h not found?

texts:
I am a Qt/C++ developer. I have to create a GUI application (Qt/C++) with GIS maps. I'm trying to do that using QGIS.

I'm working in Linux Ubuntu with Qt Creator 2.4.1 (based on Qt 4.8.0).

I did the following steps:

I installed QGIS packages (QGis starts correctly)
Then I downloaded QGIS source files.

Now I'm trying to create a simple Hello World QGIS application, but something doesn't work.

My programs has just one main.cpp file. The code is:

#include "qgsapplication.h"

int main(int argc, char ** argv)
{
   QgsApplication app(argc, argv, true);
   return 0;
}

The .pro file is:

#-------------------------------------------------
#
# Project created by QtCreator 2015-05-15T09:58:00
#
#-------------------------------------------------

QT       += core gui

TARGET = TestQGIS
TEMPLATE = app

INCLUDEPATH += /home/pesce/qgis-2.8.1/src/core/

SOURCES += main.cpp

HEADERS  +=

I'm trying to compile but the following error is shown:

/home/pesce/GUI/TestQGIS/../../qgis-2.8.1/src/core/qgsapplication.h:23: error: qgsconfig.h: No such file or directory

It seems that qgsconfig.h is missing......

Which is the problem? How can I fix this problem?

There was a libraries problem..not all the needed libraries were included!
I've included all the libraries and now all the needed files are correctly present.

But I'm not still able to compile the code.... the following error is shown:
/usr/include/qgis/qgsapplication.h:55: error: 'override' does not name a type

It seems a problem of the compiler..

Any suggestion?


texts:
I'm trying to get a message box to display some simple counts of layers and fields. Everything works fine until I add in the code to get the field count of composite layers. Below is the code for my button click event which populates the message box, and the error message I am receiving.

private void buttonCountLayers_Click(object sender, EventArgs e)
        {
            //count the point/line/poly/raster feature layers
            int groupcount = 0;
            int pointcount = 0;
            int linecount = 0;
            int polycount = 0;
            int rastercount = 0;
            string subfields = "";
            int subfieldcount = 0;

        for (int i = 0; i < layercount; i++)
        {
            ILayer pLayer = pMap.get_Layer(i);
            string name = pLayer.Name;
            pointcount = pointcount + PointLayerTest(pLayer); //returns 1 if layer has point geometry
            linecount = linecount + LineLayerTest(pLayer);
            polycount = polycount + PolygonLayerTest(pLayer);
            rastercount = rastercount + RasterLayerTest(pLayer);

            if (pLayer is ICompositeLayer)
            {
                ICompositeLayer pCompLayer = pLayer as ICompositeLayer; //QI
                int sublayern = pCompLayer.Count;
                for (int j = 0; j < sublayern; j++) //go through group layer
                {
                    ILayer pLayer2 = pCompLayer.get_Layer(j);
                    groupcount = groupcount + GroupLayerTest(pLayer2);
                    pointcount = pointcount + PointLayerTest(pLayer2);
                    linecount = linecount + LineLayerTest(pLayer2);
                    polycount = polycount + PolygonLayerTest(pLayer2);
                    rastercount = rastercount + RasterLayerTest(pLayer2);
                    IFeatureLayer2 pFeatureLayer = (IFeatureLayer2)pLayer2;
                    ITable subLayers = (ITable)pFeatureLayer;
                    subfieldcount = subLayers.Fields.FieldCount; //this line throws the error
                    subfields = subfields + Environment.NewLine + pLayer2.Name + ":" + subfieldcount;
                }
            }
        }
        //show all counts in pop-up window
        MessageBox.Show("number of group layers:"+groupcount+ Environment.NewLine +"point layers:"+pointcount+
            Environment.NewLine+"line layers:"+linecount+Environment.NewLine+"polygon layers:"+polycount+
            Environment.NewLine+"raster layers:"+rastercount+Environment.NewLine+Environment.NewLine+subfields);
    }

ERROR MESSAGE:
An unhandled exception of type 'System.Runtime.InteropServices.COMException' occurred in projectname.dll

Additional information: Error HRESULT E_FAIL has been returned from a call to a COM component.


texts: ArcGis Rest Vehicle Routing Time Problem

texts:
I am using the ArcGis vehicle routing rest service.
It works fine until I add in the EarliestStartTime for my route, then it no longer assigns the orders to the routes (delivery windows and breaks also cause the same problem).
Even if I make the EarliestStartTime 0 it still doesn't work, only removing it works.

What am I doing wrong?

request=?token=NmsVEs18OZjNkg0jAiUfX-ltrqdqXAQdoZg775KkVpPB6dPPSqm9_N0O2QgR0i8Qnz-AFbi6fGVO60_o7zD2ol43J8nKECoFlZ4QdjndfpAi5JMxtWeot7WUbD5uX-4fFggJWuOgGhiEG7zZt1hmPA..
&default_date=1356998400000
&orders={
  "features": [
    {
      "geometry": {
        "x": 115.9679,
        "y": -31.78098
      },
      "attributes": {
        "DeliveryQuantities": 1,
        "Name": "5406",
        "ServiceTime": 15
      }
    }
  ]
}
&depots={
  "features": [
    {
      "geometry": {
        "x": 115.93737,
        "y": -31.90728
      },
      "attributes": {
        "Name": "DBW593",
        "CurbApproach": 0
      }
    }
  ]
}
&routes={
  "features": [
    {
      "attributes": {
        "Name": "2",
        "StartDepotName": "DBW593",
        "EndDepotName": "DBW593",
        "EarliestStartTime": 1357002000000,
        "ArriveDepartDelay": 5,
        "Capacities": 17,
        "MaxTotalTime": 1095
      }
    }
  ]
}


texts: Arcmap - automatic print area per one feature

texts: I have hundreds of routes in vector shpfile. They have different lenghts and I want to create one map per one feature - route- in shapefile. Unfortunatelly size of print area will be different for each route (I want to keep only the same scale for each). I don't want to prepare for everyone manually size before printing. Is there tool/solution which allow to automatically prepare print size for each feature (one feature-one route) with fixed scale?


texts: Convert coordinate to tile row and tile column in Google Maps when using custom projection

texts:
I'm trying to get Google Maps (javascript SDK) to use a different tile server that serves tiles as EPSG:3006. I'm using Proj4js to convert the coordinates.

So far I have this code:

 proj4.defs("EPSG:3006", "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
 var googleProj = proj4("EPSG:4326")
 var lantmaterietProj = proj4("EPSG:3006")
 var projectionTransform = proj4(googleProj, lantmaterietProj)

 var lantmaterietMapType;
 function initLantmateriet() {
   lantmaterietMapType = new google.maps.ImageMapType({
     getTileUrl: function(coord, zoom) {
       return "https://api.lantmateriet.se/open/topowebb-ccby/v1/wmts/token/secret-token/1.0.0/topowebb/default/3006/" + zoom + "/" + coord.y + "/" + coord.x + ".png";
     },
     tileSize: new google.maps.Size(256, 256),
     minZoom: 0,
     maxZoom: 9,
     name: 'Lantmateriet'
   });

   lantmaterietMapType.projection = {
     fromLatLngToPoint: function(latLng) {
       var newPoint = projectionTransform.forward([latLng.lng(), latLng.lat()]);
       return new google.maps.Point(newPoint[0], newPoint[1]);
     },
     fromPointToLatLng: function(point, noWrap) {
       var newPoint = projectionTransform.inverse([point.x, point.y]);
       return new google.maps.LatLng(newPoint[0], newPoint[1]);
     }
   };
 }

The coord.x and coord.y in the getTileUrl function cannot be put in the url like that. I'm not sure how to calculate the correct tile row and tile column though.

Any ideas on how to calculate the correct tile row and column here?


texts: Monitor Traffic Between ArcGIS Server and Oracle DB

texts:
So, my dba informed me yesterday that ArcGIS Server made 40+ million requests to Oracle in a 20 minute period.  That would be approx. 35k requests per second.  I started laughing, and said we have the greatest performing ArcGIS Server in ArcGIS Server history - especially since its a single machine configuration.

I have no idea how he determined ArcGIS server to be making this volume of requests.  Before I go do battle with him, is there some type of quick & easy monitoring tool I can run on the ArcGIS to gauge how many requests its throwing at Oracle.


texts: How to reproject aqua modis images / ocean color?

texts:
I am working with aqua modis, ocean color, L2 - SST image but in order to do I must reproject it. I use the gdal_translate to convert the format of the image but reprojection is becoming very difficult...

When I download it, it has a north down projection. This image has 5 bands and the two first bands are the lat lon in decimal degrees.

Do you know how can I reproject them by using the band 1 and band 2 as coordinates?


texts: Automatically generating pour points along stream network?

texts:
The following screenshot shows the results of a water flow analysis where I identified stream order via the Strahler method in ArcGIS.  Now I need to generate watersheds for 2nd order streams and above. The screenshot shows the areas I need to automatically generate pour points. I have tried using the stream link tool, although the resulting product generates watersheds for all stream links, not just those >= 2nd order.

How can I automatically generate pour points for >= 2nd order streams?


texts: Finding/buying Civic Address Data for all of Canada

texts:
I work for a company that is looking to buy all the civic address data for Canada.
I've been Googling around and I have yet to find a solution that is not pay per request.

We need it so that it can be integrated in to an ArcGIS server.

What we want to do with the data is tie to in to a Geodatabase where people will be able to search their address and have information pop up about their house. We have the metadata for each location but our current addressing system is partially interpolated/inaccurate (not accurate to the level that we'd like).

We're looking for source data.


texts: Showing only selected country from GeoJSON file and hiding remaining unselected country using OpenLayers?

texts:
How do I not show any countries except Zambia?

<script type="text/javascript">
         (function() {
         var inputs  = 'Zambia';//may be any country
         var raster = new ol.layer.Tile({
            source: new ol.source.OSM(),
           projection: 'EPSG:3857'
          });
        var vector = new ol.layer.Vector({
             source: new ol.source.Vector({
             ormat: new ol.format.GeoJSON(),
              url: 'countries.geojson'
            }),
           style: function(feature, res){
               if(feature.get("name") == inputs){
                  return new ol.style.Style({
                  stroke: new ol.style.Stroke({
                   color: 'red',
                   width: 2
                  })
                   });
              }
          }
       });
     var map1 = new ol.Map({
        target: 'map',
        layers: [raster, vector],
        view: new ol.View({
     center: ol.proj.transform([27.60315,-13.41099],'EPSG:4326','EPSG:3857'),
     zoom: 6
    })
     });
  })();
</script>


texts: Rotating Layout View and Labels in ArcMap?

texts:
I have a question if you are able to, upon rotating the map, is there a way to do the same with the labels? If so, how?

Also, besides using the "Rotate Data Frame" to rotate the Layout View, which other method can you do to rotate the layout? I right-click but the rotate options are ghosted, meaning those options are non-clickable.

I have ArcView 9.3.


texts:
Using ArcGIS 10.2.2
I am trying to set the extent in environment settings using arcpy.env.extent

My code is here:

for lyrVec in arcpy.mapping.ListLayers(mxd, "", df):
    if lyrVec.name == "Thresholds":
        desc = arcpy.Describe(lyrVec)
        lyrVecPath = desc.path + '\\' + str(lyrVec) + '.shp'
        fixedpath = os.path.normpath(lyrVecPath)
        arcpy.env.extent = lyrVecPath

Error message is thrown up:

  RuntimeError: Object: Error in accessing environment


texts: Unable to select SQL Spatial table line feature

texts:
I have a table in a SQL database that has been spatially enabled with a geography column.

In ArcGIS it is a line feature and can be viewed in ArcMap with no issues and even published as a map service with ArcGIS for Server.

However, neither the feature when it is in ArcMap or as a map service in a web app is able to be selected through the map. It can be selected through the attribute table.

I have adjusted and tested all selection settings in ArcMap. I have exported the data out as a shapefile an it works correctly. I have another table in the same SQL database that is a point feature and it works correctly.

Has anyone experienced this before or know of any esri bug related to this?

UPDATE 1:

I have rebuilt the table using geometry instead of geography - no change.

There is one feature that can be selected, it is the only feature west of the prime meridian.

UPDATE 2 (in response to comments):

1- Not sure about spatial indexes with SQL Spatial tables, please advise what i should be looking for?

2- I added the table to ArcMap from the database connection made in ArcCatalog, it calculated the extent on import, and maintained the correct spatial reference information as below:

  GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

The data also draws in the correct locations.

3- The table does have an objectID field. The table is not in a geodatabase, it is in SQL Server.

UPDATE 3:

Have found the below known issue. Not exactly what I am experiencing but similar. I adjusted the data to only be inside one hemisphere and it is still not working.

Bug: Spatial query fails...

UPDATE 4:

Have tested the same data in Blue Marble Global Mapper, and it works correctly. Looking like an esri bug.


texts: Convert shapefile to KML without losing symbols

texts:
I am using QGIS 2.2 and want to convert a shapefile containing point coordinates with userdefined svg-symbols to a KML file.

Everything seems to work. The original file is in WGS 84 / UTM zone 32N and I convert it to KML EPSG:4326. When I open the new KML file in Google Earth the position of the points is correct but as a symbol I get Google Earth-Pins.

Is it possible to include the svg-symbol and the label as well in the KML file?


texts: Leaflet - Load markers in the current map - dynamically

texts:
I am looking for a marker plugin for leaflet that allows loading of markers dynamically, according to the current view.

E.g. I want to add OSM POIs to a world-wide map. Obviously loading all at once would be too much. Therefore, I'd only load markers on smaller areas when zoomed in. Is there a plugin that supports this functionality?

I know there is one for the overpass-api (https://github.com/kartenkarsten/leaflet-layer-overpass), but I do not want to spam the Overpass Api. Also the response is to slow for a customer facing application (IMHO). Additionally, I do not want to only map OSM POIs.

I thought about setting up my own DataBase and changing the leaflet-layer-overpass plugin to work with that, but that seems to be a lot of work. Is there maybe a pre-existing solution to the problem?


texts: What is the difference between the compact build and standard build of the ArcGIS JavaScript API?

texts:
I am creating a JSAPI application using version 3.2/3.3 and have recently gone through the page detailing the difference  between the compact build and standard build.

Right now, I am confused about which build I should use. (I know a few people will use iPhone while others will use a PC to access our application.) So what should I use in my application?
Since in the compact build, a few tools are not supported, I cannot create all applications in it.

Software platform: ArcGIS Server & Desktop 10.1

Device platform: Windows , iOS & Android


texts: someone know how to abilitate ECW support on QGIS version 2.14.4-Essen (Ubuntu 14.04)?


texts: Correct folders and conf settings for qgis-web-client with qgis-server?

texts:
I have setup QGIS server using this guide sucessfully and would now like to setup qgis-web-client alongside it.

My QGIS server folder is: /usr/lib/cgi-bin

For qgis-web-client I have followed the qgis-web-client readme here and have cloned from github (as recommended) here:

/home/web/qgis-web-client

I am now stuck with editing the qgis-web-client.conf file to get things working correctly. I updated the path to "home/web/qgis-web-client/site" file. (See qgis-web-client.conf below for reference)

I think I am doing something wrong now and cannot understand the rest of the instructions in the readme. In particular:

  Check the symlink in cgi-bin is correct.
  Check the QGIS libs are in your /etc/ld.so.conf path
  Copy site/index.xml and check paths match your system OR
  Modify index.html and point your browser to that

Please can you tell me what is wrong with my configuration?

Please can you explain what the 3 points above mean?

Many thanks!

There seems to be many questions related to how to setup qgis-web-client and a lot of confusion. I have spent a lot of time scouring the gis stackexchange for the answer but so far have been unsuccessful!

<VirtualHost *:80>
ServerAdmin webmaster@localhost
# Add an entry to your /etc/hosts file for qgis-web-client.localhost e.g.
# 127.0.0.1 qgis-web-client.localhost
ServerName qgis-web-client.localhost

DocumentRoot /home/web/qgis-web-client/site
<Directory />
Options FollowSymLinks
AllowOverride None
</Directory>
<Directory /home/web/qgis-web-client/site/>
Options Indexes FollowSymLinks MultiViews
AllowOverride None
Order allow,deny
allow from all
</Directory>

ScriptAlias /cgi-bin/ /home/web/qgis-web-client/cgi-bin/
<Directory "/home/web/qgis-web-client/cgi-bin/">
AllowOverride None
Options +ExecCGI -MultiViews -SymLinksIfOwnerMatch
Order allow,deny
Allow from all
#Uncomment the next line to enable logging to a file
#SetEnv QGIS_LOG_FILE /tmp/qgislog.txt
</Directory>

# optional rewrite rules

# Forbid direct access
#  RewriteRule ^/cgi-bin/.*$ - [F]

# Search with SearchPanel

# sample search on layer 'Hello'
RewriteCond %{QUERY_STRING} ^(?:.*)query=samplesearch&*(?:.*)$
RewriteCond %{QUERY_STRING} ^(?:(?:.*)&)?colour=([^&]*)(?:.*)$
RewriteRule ^/wms/(.+)$ /cgi-bin/qgis_mapserv.fcgi?map=/<path-to-qgis-server-      projects>/$1.qgs&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&LAYERS=Hello&QUERY_LAYERS=Hello&FEATURE_COUNT=20&INFO_FORMAT=text/xml&SRS=EPSG:4326&FILTER=Hello:"colour"\ =\ '%1' [PT]

# Rewrite /wms/mapname to qgis_mapserv.fcgi?map=/<path-to-qgis-server-    projects>/mapname.qgs
RewriteRule ^/wms/(.+)$ /cgi-bin/qgis_mapserv.fcgi?map=/<path-to-qgis-server-projects>/$1.qgs [QSA,PT]
# Rewrite /maps/mapname to qgis-web-client main page. mapname will be extracted for wms calls in Javascript code.
RewriteRule ^/maps/([^\.]+)$ /qgis-web-client/site/qgiswebclient.html [PT]
# Rewrite /maps/* to qgis-web-client/site (e.g. /maps/gis_icons/mActionZoomNext.png -> /qgis-web-client/site/gis_icons/mActionZoomNext.png)
RewriteRule ^/maps/(.*) /qgis-web-client/site/$1 [PT]

ErrorLog ${APACHE_LOG_DIR}/qgis-web-client-error.log
CustomLog ${APACHE_LOG_DIR}/qgis-web-client-access.log combined

</VirtualHost>


texts: Raster is flipped when opened as numpy array

texts:
I have two rasters of the same dimensions. Opened in QGIS, they align perfectly. However when I open them via GDAL as a numpy array, one of them is flipped and I have to flip it back to make calculations on it. See my code:

ds = gdal.Open('{}/data/tmp/one25832.tif'.format(ortascii))
band = ds.GetRasterBand(1)
arrA = band.ReadAsArray()
arrA=np.flipud(arrA)

file1 = '{}/data/tmp/two25832.tif'.format(ortascii)
ds1 = gdal.Open(file1)
band1 = ds1.GetRasterBand(1)
arrB = band1.ReadAsArray()

arrB = np.where((arrA ==999999), arrB, arrA)

Any insights on how this can happen? They have the same CRS...


texts: Size scale field - QGIs Pisa 2.10-1 Where did it go?

texts: I've been scaling symbol sizes in QGIS using the size-scale field under the "advanced" option on the properties-style dialog.  We just upgraded to Pisa and now it's gone - anyone know where it went?  I'm tempted to roll back to 2.8 if I can't find it.


texts: Selecting Maximum value for each unique ID

texts:
So i have a table (linked below) and for each unique PropertyID, i need to find the highest percentage value and extract that to a different table. I need to retain all the other attribute information as well, (which summarize or dissolve doesn't do).

I'm using ArcMap 10.3


texts: leaflet-> draw a rectangle with center point

texts:
In leaflet, I would like to draw the largest sqare that fits inside a circle. In other words, I have the radius info and the center, but not the corner points.
How can I achive that?

And here is my code for that circle (cyc_163_1):

var clausis_163_cyclic =  L.circle([44.6730555555556, 6.9225], 700, {fillColor:'none', color: '#a6a6a6', weight: 6, opacity: 0.3});
var length = 0.0065
var angle = 254    //16juillet->270-16=254
var end_x = clausis_163_cyclic.getLatLng().lng + length * Math.cos(angle * Math.PI / 180)
var end_y = clausis_163_cyclic.getLatLng().lat + length * Math.sin(angle * Math.PI / 180)
var cyc_163_1 = L.circle([end_y, end_x], 80, {fillColor:'#666666', fillOpacity: 1, color: '#666666', weight: 0.5, opacity: 1});

Normally I believe that the following code should solve the issue, but it does not:

var rectangle = L.rectangle(circle.getBounds(), {
  stroke: false,
  fill: true,
  color: '#00f',
  opacity: 0.5
}).addTo(map);


texts: No Layers available in GeoExt Styler for GeoServer?

texts:
I´m searching for a styler to create styles for my data and I stumbled over the styler extension for GeoServer 2.6.1.

I that the REST-Plugin is now is integrated in GeoServer and i only had to "install" the styler-plugin (extract the zip-file to the www-directory).
If I open the styler-page (localhost:8080/geoserver/www/styler/index.html) in the browser I get the styler-page but the layer-panel is blank. There are no layers shown.
It seems that there is no connection between the styler and the geoserver.

Is there anything i have to configure somewhere?


texts: GeoDjango Setup: Python 'GDALAllRegister' Attribute Error

texts:
I am trying to setup Geodjango 1.8 in a virtualenv with Python 2.7.9 (Win32).

I have installed and setup Postgres/PostGis and have my Path setup to my OSGeo4W directory. I have also installed the required Python bindings.

When I try to 'runserver' (to start the development server), I get the follow error:

    ...
    from django.contrib.gis.gdal.prototypes import ds as vcapi, raster as rcapi
    File "C:\Users\Matt\Envs\django_venv\lib\site-packages\django\contrib\gis\gdal
    \prototypes\raster.py", line 24, in <module>
    register_all = void_output(lgdal.GDALAllRegister, [])
    File "C:\Python27\Lib\ctypes\__init__.py", line 378, in __getattr__
    func = self.__getitem__(name)
    File "C:\Python27\Lib\ctypes\__init__.py", line 383, in __getitem__
    func = self._FuncPtr((name_or_ordinal, self))
    AttributeError: function 'GDALAllRegister' not found

I am able to start the development server in an app that doesn't use 'django.contrib.gis'.

My Settings.py:

DEBUG = True

ALLOWED_HOSTS = []

# Application definition

INSTALLED_APPS = (
'django.contrib.admin',
'django.contrib.auth',
'django.contrib.contenttypes',
'django.contrib.sessions',
'django.contrib.messages',
'django.contrib.staticfiles',
'django.contrib.gis',
'leaflet_map',
'djgeojson',
)

MIDDLEWARE_CLASSES = (
    'django.contrib.sessions.middleware.SessionMiddleware',
    'django.middleware.common.CommonMiddleware',
    'django.middleware.csrf.CsrfViewMiddleware',
    'django.contrib.auth.middleware.AuthenticationMiddleware',
    'django.contrib.auth.middleware.SessionAuthenticationMiddleware',
    'django.contrib.messages.middleware.MessageMiddleware',
    'django.middleware.clickjacking.XFrameOptionsMiddleware',
    'django.middleware.security.SecurityMiddleware',
)

ROOT_URLCONF = 'webmaps.urls'

TEMPLATES = [
    {
        'BACKEND': 'django.template.backends.django.DjangoTemplates',
        'DIRS': [],
        'APP_DIRS': True,
        'OPTIONS': {
            'context_processors': [
                'django.template.context_processors.debug',
                'django.template.context_processors.request',
                'django.contrib.auth.context_processors.auth',
                'django.contrib.messages.context_processors.messages',
            ],
        },
    },
]

WSGI_APPLICATION = 'webmaps.wsgi.application'

# Database
# https://docs.djangoproject.com/en/1.8/ref/settings/#databases

DATABASES = {
    'default': {
        'ENGINE': 'django.contrib.gis.db.backends.postgis',
        'NAME': 'db',
        'USER': '',
        'PASSWORD': '',
        'HOST': '',
        'PORT': '',
    }
}

# Internationalization
# https://docs.djangoproject.com/en/1.8/topics/i18n/

LANGUAGE_CODE = 'en-us'

TIME_ZONE = 'UTC'

USE_I18N = True

USE_L10N = True

USE_TZ = True

# Static files (CSS, JavaScript, Images)
# https://docs.djangoproject.com/en/1.8/howto/static-files/

STATIC_URL = '/static/'


texts: Splitting multilinestrings in segments using PosgreSQL

texts:
I need to split several multilinestring into a defined number of segments using PosgreSQL.

I tried using ST_subdvideand ST_segmentize but they didn't work.

Can someone help me to solve this problem?

I tried using the following SQL string to split the first line in segments of 2 m:

with x as(

select * from city.street t
where id_street = 1)

select
id_street,
st_astext(geom),
st_astext(st_segmentize (x.geom, 2))from x;


texts: How to export Selection Feature to CAD in ArcGIS 10.2.2

texts:
This code can not be exported select feature to CAD file

ESRI.ArcGIS.ConversionTools.ExportCAD exportCad =  new ESRI.ArcGIS.ConversionTools.ExportCAD();

exportCad.in_features = @"D:\mydb.gdb\Street";

exportCad.Output_File = @"C:\StreetCAD.DWG";

exportCad.Output_Type = "DWG_R2007";

exportCad.Append_To_Existing = "true";

Geoprocessor gp = new Geoprocessor();

gp.Execute(exportCad, null);

I want to export only feature selection, who can help me!


texts: How to detect overlapping linestrings in PostGIS?

texts:
I'm looking for the fastest way to detect whether a given linestring overlaps another linestring. (As reference, let's say my wife and I both go for a run and take different routes, I'd like to detect the segments our different runs had in common.) Here's my current approach.

Given two tables: r1 and s1, both containing a geo_data column of linestring, first narrow the data (taking advantage of gist index) with r1.geo_data && s1.geo_data, then use ST_Contains to see if s1 is overlapping r1, i.e.

r1.geo_data && s1.geo_data and ST_Contains(ST_Buffer(r1.geo_data,ST_Length(r1.geo_data)*.002),s1.geo_data)

This approach works, but it's extremely CPU-intensive and slow. Is there a better approach?


texts:
I have 2 rasters, one representing land-cover and 1 representing trails and roads. I want to use the con function in ArcMap so that my new raster will contain both but the unclassified or wrongly classified pixels on the landcover for roads will be replaced due to the trails/roads raster layer.

I have read the documentation but dont fully understand which raster layer is the input conditional and which is the input true raster? Or how to construct the query? I know what i want to say in english but I dont know how to write it in SQL.


texts: Accessing the features of a temporary output with PyQGIS

texts:
I want to access the features of a temporary output to save certain geometric properties (in this case the length) in a variable.

Usually I would use:

path="C:\Users\..."
for files in os.listdir(path):
    if files.endswith(".shp"):
        layer = QgsVectorLayer(path+'/'+files, files, "ogr")
        features = layer.getFeatures()
        for f in features:
            geom = f.geometry()
            length = geom.length()

But since I'm now working with a temporary layer, created in a previous step in the script, I want to do something like this:

path="C:\Users\..."
for files in os.listdir(path):
    if files.endswith(".shp"):
        output=processing.runalg("qgis:pointstopath",path+'/'+files,"field1","field2","",None,None)
        features = output['OUTPUT_LINE'].getFeatures()
        for f in features:
            geom = f.geometry()
            length = geom.length()

How can I access the features of the temporary output['OUTPUT_LINE'] layer?


texts: How to Get Dynamic Map Service Layer Name On ArcGIS API

texts:
Can you please take a look at This Demo and let me know how I can get the Selected Layer Name through API when check boxes are checked?

As you can see at This Link the service is providing 3 layers as:

and I need to get these names as they are selected
I already tried the Request layer info but this seems retuns metadata about fields only.