texts: How to create cached service from 2TB of .TIFFs to ArcGIS Server without buying Imagery Extension?

texts:
I work for a local government enterprise GIS team in California. We are an ESRI shop and were able to piggy back on our county's grant-funded aerial imagery purchase this year. We didn't have to go through the normal procurement process and were essentially gifted a bunch of high resolution imagery for our city.
I've been having issues coming up with a least cost way to get these published on our ArcGIS Server for sharing with others.

I've got 3 different formats of imagery: TIFF, JP2 and MrSID
The JP2s look horrible so we essentially have 2 formats. Our initial plan was to use ArcGIS to Mosaic the TIFFs and then create a service to be cached upon publishing. After going through the process, we found out that we'll have to purchase the Imagery Extension to be allowed to publish this way. Imagery Extension is ridiculously expensive, so I figured I'd see if I could find alternatives.

Has anyone had success in working around the Imagery Extension to cache a large amount (2TB) of imagery? I've found examples of ways around this but nothing on this scale. ArcMap crashes if I add half the images to the map.


texts: QGIS - export a project with dependencies

texts:
I was wandering if there was a way to save a specific project with all layers as dependent files to make it easy to share with collaborators.

for example, create a zipped or non-zipped folder containing the same structure as the project :

Root_folder/
my_project.qgs
    \____________myfolder
                     \___Layer 1
                     \___Layer 2
                     \___Layer 3
    \____________myfolder2
                     \___Layer 4
                     \___Layer 5
                     \___Layer 6
    \____________myfolder3
                     \________myfolder4
                                     \___Layer 6
                                     \___Layer 7
                     \___Layer 8

If it does not exists, would it be tricky to develop regarding the current API ?


texts: Which interface controls 'Place overlapping labels' for a feature layer?

texts:
I need to turn on the option to 'Place overlapping labels' for a feature layer via ArcObjects. In the ArcMap UI it is located here:

I can get to most of the labeling options from the IAnnotateLayerProperties, ILabelEngineLayerProperties and IBasicOverposterLayerProperties interfaces.

The option to control overlapping labels seems to be missing.


texts: It is possible to use different fields for symbology and legend labeling in QGIS?

texts:
[Edit : working with Quantum GIS version 2.0.1-Dufour]

I wonder if it is possible to use some field to set the symbology and some other to label the map.

[Edit 1 : I mean labeling the legend in the map designer tool to output as JPG or PDF, not labeling objects in the canvas for working purpose.]

[Edit 2 : The purpose of my question is to get an automatical process so I could save time and not manually change all the labels of the legend.]

My problem is : I have to make a map of natural habitats, based on european "CORINE biotopes" typology which uses entities with number and names, e.g. :

"37 hygrophilous meadows",
"38 mesophilous meadows",
"41 deciduous forests",
"42 coniferous forests",
"51 active bogs",
etc.

So I would like to use shemes of the same color for habitats which are similar : all meadows in different tons of yellow, all forest in different tons of green, mires and bogs in some brown scale etc.

To do this I need to base my styling work on the NUMBERS, since following numbers identify similar habitats (meadows are 3X, forests are 4X, etc).

But for reading and understanding purpose, the printed map should be labelled with the NAMES. And the labelling should not get "noised" by the numbers, so I separated NAMES and NUMBER of the same entities in 2 different fields.

Is there some way to say QGIS to let me do the symbology on the basis of the "NUMBER"-field, but then to label the map on the basis of the "NAME"-field ?

I would like to set the "label" (like here from the Layer > Properties > Symbology)

To set the label of the legend in the Map Designer Tool : Here I would like to get the names of the habitats, and not the numbers.

Thanks for all kind of help.


texts: How to translate GDAL Terminal code for python use

texts:
I am trying to align two raster image with each other using Gdalwarp.The code below  is Gdalwarp code.
gdalwarp  -tap -t_srs epsg:32614 F:\\1.tif F:\\out1.tif
gdalwarp  -tap -t_srs epsg:32614 F:\\2.tif F:\\out2.tif

I translate this code for python use
from osgeo import gdal

input_file = gdal.Open('F:\\2.tif')
ds = gdal.Warp('F:\\output2.tif', input_file, dstSRS='EPSG:32614',targetAlignedPixels= True)
ds = None

This code runs but is not giving me output


texts:
I am hoping that someone can point me in the right direction. I am wanting to calculate the distances from a point to each point in another layer along a road network using v.net.distance. I have done this successfully in the past. However I now cannot get it right. All layers are in the same projection.

The steps I have taken:

Create starting point layer
Create destinations point layer
Clip road network to area of interest
Use v.net to connect starting point to road network
Use v.net to connect destination point to road network
Use v.clean (snap and break) to clean topology
Use v.net.distance to calculate shortest path distances

I am now getting an error that the layer is not correctly generated:
Loading resulting layers

  "The following layers were not correctly generated. .../dist.shp You can check the 'Log Messages Panel' in QGIS main window to find more
  information about the execution of the algorithm."

Inspection of the log shows the following error:

  "ERROR: Value  out of range for parameter  Legal range: point,line,boundary,centroid,area,face,kernel,auto "


texts:
This question was asked before on Stack Overflow.

Let's say we have mercator map with given resolution (for simplicity in centimeters).

Now let's get coordinates of (lat,lon) = (50, 15):

$ proj +proj=merc +ellps=WGS84
50 15
5565974.54      1678147.52

Since we don't know the boundaries, we can't scale it to desired resolution.

The only way I've found to get boundaries is compute the coordinates for (lat, lon) = (180, 0) and (0, 90), but that's more like a hack and it's not universal (other projections might map only from (±90, ±90)).

Does proj4 deal with scaling differently?


texts: Getting geo transforms (affine) without a TIFF reference image

texts:
I have a binary mask with buildings' footprints.
What I have is a latitude and longitude of the center of my image and I want to georeference each polygon first.
And here is what I have done before:

I get the contours (pixel coordinates) for each polygon and store them in a data frame using OpenCV
 contours = measure.find_contours(self.bin_mask, 0.5).

Now I need to get the georeferenced coordinate. I have already done that before with rasterio:
(px, py) = rio.transform.xy(self.tif_img.transform, x, y, offset='center')

Where I get the georeferenced coordinates for each pixel in a polygon
Everything was working great and I was able to extract the building, the issue now is that I don't have the tiff image.
Is it possible to get the geo affine (transform) without the referenced TIFF?
In short, I need a way to georeferenced the image pixels and all I have is:
1- a lat/long coordinates of the center of the image.
2- The images size 650X650 and 1 pixel = 0.297 m
3- coordinates is (-95.407189656922 29.893852182227928)
Update
1- I know that I need to get the bounding box for my image. From what I understood so far, this will help me get the geo transform
2- I also tried this solution but it didn't work, I think I need more than one GCP to implement this correctly.


texts: ArcGIS crashes when trying category symbology?

texts:
I'm facing an issue with ArcGIS 10.3. When I try to use symbology type "category" the program crashes. It seems that is not a matter of data because I used the data before and nothing happened. The problem started yesterday after a windows 10 update so I consider something relate do this update (may 10).

I tried the complete uninstall and reinstall and didn't solve the problem.

I tried also save the MXD with another name in other location.

So my question is related to the possibility of a problem caused by an windows update.

Any ideas?


texts:
I have this tool where the user inputs an ID number and the tool selects a list of features applies symbology (.lyr) and adds the queried features to the open mxd.

# Import modules
import arcpy
import os

# Set overwrite option
from arcpy import env
env.overwriteOutput = True

# Set variables
pws = arcpy.GetParameterAsText(0)
files = ["SERVICE_AREA_COMBO", "BACTI", "DBP", "DBP_STAGE2", "OFFICE", "PRV", "PUMP_STATION", "TANK", "VALVE_MISCELLANEOUS", "WTP"]  # list all names used for feature classes and .lyr files.
layers = ["Service Area " + pws, "Bacti " + pws, "DBP " + pws, "DBP Stage 2 " + pws, "Office " + pws, "PRV " + pws, "Pump Station " + pws, "Tank " + pws, "Valve " + pws, "Water Treatment Plant " + pws]
root = r"V:\WATER_FACILITY\WATER_FACILITY.gdb"
symRoot = r"V:\ArcGIS_Tools\pws_id_tool"

for i, file in enumerate(files):
    path = os.path.join(root, files[i])  # Create directory path for featureclasses
    symPath = os.path.join(symRoot, files[i] + ".lyr")  # Create directory path for .lyr files
    arcpy.MakeFeatureLayer_management(path, layers[i], "PWS_ID" + " = '" + pws + "'")  # Make all the layers from "files" list
    arcpy.ApplySymbologyFromLayer_management(layers[i], symPath)  # Apply the symbologies from "symNames" list.
    featureCount = int(arcpy.GetCount_management(layers[i]).getOutput(0))
    arcpy.AddMessage(layers[i] + " count = " + str(featureCount))
    if featureCount > 0:
        newLayer = arcpy.mapping.Layer(layers[i])
        arcpy.mapping.AddLayer(df, newLayer, "TOP")
    else:
        arcpy.AddWarning("No " + layers[i] + " features for this PWS.")

The only parameter entered is ID number.  I would like to have the user select which features they want to add to the map by using GetPramaterAsText  with the 'Data Type' set as boolean.  That way they can check, or uncheck, the boxes based on what features they want added to their mxd.  I' m not quite sure how I would approach this with the for loop, if it can be done?


texts: Changing the scale (generalisation) in ArcGIS 10

texts:
I have 2 polygon layers, accordingly 1:10 000 and 1:500 000.
I want to generalize the first layer from 1:10 000 to 1:500 000.

I tried using Simplify polygon tool, and since these scales relate to each other as 1:50, I used these options in the tool:

Simplification Algorithm : POINT_REMOVE
Maximum allowable Offset : 50 Meters

Is this the right approach, because the result does not look like what I was looking for. There is only a subtle change in detail level.


texts: Getting QgsMapCanvas viewport's four corners' coordinates

texts:
I have a PyQGIS application that uses the QgsMapCanvas as map display. I am trying to retrieve its four corners' coordinates, but couldn't find any public function to do that. All I have is its width() and height().

Is there any way to get the actual coordinates of the corners of the viewport?

Let me provide further context, I need the map canvas's viewport's corners' coordinates when it has not had any layers loaded into it yet or the first layer I am constructing is in memory that has yet to set its projection, I found that in such cases, the map canvas's extent is empty, so calling its extent won't help.


texts: Can 3rd party javascript libraries be used with CartoDB?

texts:
I'd like to build some functionality (search, sliders, etc.) in a 3rd party javascript library (jquery mobile), and integrate it with CartoDB.

Can this be done?  Are there any examples available or a cookbook decribing how to do the integration ?


texts: Convert CensusWard codes to Google Maps (UK)

texts:
The UK's Office for National Statistics uses Census Ward codes for some of its data.

For instance, this file - https://www.ons.gov.uk/visualisations/dvc479/map/data/data.xls - has life expectancy by ward.

Example codes:

E36000489 Bathavon North
E36000492 Bathwick

What's the best way to turn these codes into a map showing the border of the area?


texts: Calculating rate of spread of tree species

texts:
I'm fairly new to GIS and using ArcGIS for a palaeoecology project at university. I have created a map with isolines (all in the same layer using the contour tool) which represent time of arrival of a tree species.

I need to work out the rate of spread in km per year but have no idea how to do this!
I presume I have to find out the average distance travelled and divide by the time taken but I'm not sure how to measure the distance, as it has spread in more than one direction.  Does anyone have any ideas? I've spent the last two weeks trying different tools and getting nowhere!


texts: Using TLE data how would I create a KML file to look at a particular spot at a particular time

texts:
Let's say I have some TLE data as follows:

REDACT
1 13337U 98067A   20087.38052801 -.00000452  00000-0  00000+0 0  9995
2 13337  51.6460  33.2488 0005270  61.9928  83.3154 15.48919755219337

I wish to view, for example, the Washington monument on 1st June 2019 14:00

I'd want this in a KML file so I can import into Google Earth Pro.

Any advice?

CLEARING UP

This is a competition at work over this holiday weekend
I have a KML file and I essentially need to fill in the blanks under NetworkLink > Link
I have the TLE information above
An image was captured from that satellite that shows a picture at a specific lat/lng at a specific time
I presume I have to see what position that satellite was at that specific time and then plug in the relevant files into the KML file (as below)

<viewFormat>BBOX=[bboxWest],[bboxSouth],[bboxEast],[bboxNorth];CAMERA=[Long],[Lat],[lookatRange],[lookatTilt],[lookatHeading];VIEW=[horizFov],[vertFov],[horizPixels],[vertPixels],[terrainEnabled]</viewFormat>


texts: How to keep the same name in the entry of the script (QGIS 2.14) output file?

texts:
For example, I have a raster called "raster_img.tiff". I use it as input in QGIS script and want it to be in output file called the same name like input + new.

##Input_raster1=raster
##Band_raster1=output raster
from qgis.core import *
from qgis.gui import *
from qgis.utils import iface
from PyQt4.QtCore import *
from PyQt4.QtGui import *
import os
from osgeo.gdalconst import *
import numpy
import sys
import osr
driver = gdal.GetDriverByName("GTiff")
ds = gdal.Open(Input_raster1, GA_ReadOnly)
dst_ds_red = driver.Create(Band_raster1, ds.RasterXSize, ds.RasterYSize, 1, GDT_Float32)
dst_ds_red.SetProjection( ds.GetProjection())
geotransform = ds.GetGeoTransform()
if not geotransform is None:
  st_ds_red.SetGeoTransform( geotransform )

bandList = []
band = ds.GetRasterBand(2)
data = band.ReadAsArray()
bandList.append(data)

rasterlayer_red = numpy.zeros((ds.RasterYSize, ds.RasterXSize), dtype=numpy.float32)
for i in range (0, ds.RasterXSize):
  for j in range (0, ds.RasterYSize):

    RED = (float(bandList[0][j][i]))
    rasterlayer_red[j][i] = RED
dst_ds_red.GetRasterBand(1).WriteArray(rasterlayer_red)


texts: Measure area of a polygon in specific units using ArcPy

texts:
I am trying to find the area of a particular polygon (OriginalPolygon, below), and I can do that... But, I'd like it in acres regardless of the original units of the polygon.(Data is coming into this script in a range of projections, but we're always interested in "acreage.")

If I was calculating geometry to store in a field, then I could use feat.area@acres -- but I don't need this stored as an attribute, I just need it as a temporary variable to tell the user. This question implies that unless I'm actually calculating geometry, the units of the shape area measurement are tied to the projection and therefore fixed.

So: Is there any way to calculate geometry in such a way that the value is only a numerical output that I can display as the script runs, or am I stuck with making an attribute field, calculating geometry, reading the field value to tell the user, and then deleting the field?

rows = arcpy.SearchCursor(OriginalPolygon)
shapeName = arcpy.Describe(OriginalPolygon).shapeFieldName
for row in rows:
    feat = row.getValue(shapeName)
    PolygonArea = feat.area


texts: postGIS SQL request with aggregating ST_Union()

texts:
I have a postGIS table that looks like this :

id  |  value  | geometry(MultiPolygon)
1      10       ...
1      20       ...
1      30       ...
2      10       ...
2      20       ...
2      30       ...

What I wish to do is to create a new geometry column containing a Multipolygon which is the union of all the other polygons of rows having the same id and a smaller or equal 'value'.

So :

first row 'new_geometry' would be identical to 'geometry'.
Second row 'new_geometry' would be the union of the first and second row 'geometry'.
Third row : union of rows 1,2,3
equal to row 4
etc

I have written something like this:

SELECT t.id, t.val, St_Union(t.geometry)
FROM
    (SELECT cs.id, cs.value as val, cs.geometry
     FROM
         mytable cs
         INNER JOIN
             (SELECT DISTINCT value from myTable) v
             ON (v.value <= cs.value)
     ) AS t
GROUP BY
    t.val, t.id;

But the geometry column I get is not aggregated, it is equal to the original geometry column.
The aggregation work however if I 'GROUP BY' only by one criteria.


texts: How to make two raster images with different column and row values the same?

texts:
I need to use two raster maps that contain exactly the same number of rows and columns to generate a relationship (For example scatter plot of Evapotranspiration (ET) and EWTcanpoy in my study area. For now, two raster images do not have the same number of rows and columns.
For example, I have 211 columns and 142 rows in the ET raster image and 1009 columns and 490 rows in the EWTcanopy raster image.
However, I adjusted my cell values same for both layers (x,y with 0.00416666667) each and it is the same for both layers. But I have different columns and rows in rasters.
How can I make both in to a same level.
Please see the attached doc for the raster information that I have.


texts: How to collect tweets outside a given bounding box?

texts:
I checked the twitter API here.

I know how to get tweets in United States within a bounding box such as:

[-124.47,24.0,-66.56,49.3843].

How can I easily get all the tweets NOT in this bounding box?


texts:
I have a number of layers which share the same attribution. When more then one of them are visible at the same time I get something like this in my attribution control:

  Company1 hyperlink1 Company2 hyperlink2 Company1 hyperlink1
  Company1 hyperlink1

How can I avoid duplicates and have only unique records shown?


texts: Adding selected feature to new layer using PyQGIS?

texts:
I'm trying to isolate one or more features in a layer, and add them to a different one. I've set up an empty layer (tempLayer) with the same fields as the one I'm selecting features from (layer1).

I then use:

query = '"AreaIDlvl0" = 25007'
selection = layer1.getFeatures(QgsFeatureRequest().setFilterExpression(query))
layer1.selectByIds([k.id() for k in selection])
for feature in layer1.getSelectedFeatures():
    tempLayer.addFeatures([feature])

But it returns False.
I've also tried adding all features:

for feature in layer1.getFeatures():
    tempLayer.addFeatures([feature])

it returns False for all of them.


texts: Map gridded rainfall data to lat-longs and extract daily rainfall data?

texts:
My objective is to create a dataset, wherein for certain locations (lat-longs) I need to report daily rainfall from 1950 to 2017.

Rainfall data source: Indian Meteorological Department (IMD) daily rainfall file - in .grd format from 1950 to 2017.

Input file: a .csv file with Latitude and Longitude columns for each location of interest.

Logic: I need to map the lat-longs to a grid (size=25km^2) and extract the rainfall on each day i.e. 01/01/1950 to 31/12/1950 from a single 1950 rainfall file. Then repeat this exercise for all the years, for the same coordinates.

I need help for starters understanding how I can see rainfall data at a daily level in a dataframe of sorts i.e. deal with .GRD files with gridded rainfall data. The software I'm hoping to use is R, because my peers have previously mapped the same rainfall data to shape-files i.e. polygons of area (district/village etc.)

I'm including the code they used for reference, this is what I was trying to use but with an input file of lat-longs.

Reference code:

*comment:Load libraries

library(dplyr)
library(rgdal)
library(sp)
library(raster)
library(lubridate)
library(magrittr)
library(readr)
library(bigmemory)

*comment:load shape file
shp<-readOGR(".","shapefile") # .shp extention not needed
myconn <- file("ind2017_rfp25.grd","rb") # read any binary file
a <- readBin(myconn, numeric(), n=135*129, size=4, endian = "little")

close(myconn)
##### start template creation #######
outfile<-paste("p2017_1","txt",sep=".")
cat(paste("ncols         135",
          "nrows         129",
          "xllcorner     66.5",
          "yllcorner     6.5",
          "cellsize      0.25",
          "NODATA_value  -999 \n",sep="\n"),file=outfile)
cat(a,file=outfile, append=T)

r<-raster("p2017_1.txt") ## reads the newly created ascii file as a raster
r<-flip(r,direction='y')
proj4string(r)<-proj4string(shp) ## assigns a coordinate system
##### end template creation #######

*comment:set no. of columns and rows according to resolution (info should be there in readme file)
n.cols <- 135
n.rows <- 129

binfiles <- dir(pattern = "ind")

*comment:1951 to 2017
*comment:leap year = 17 , non-leap = 49

l.years <- 8 # no. of leap years in the data
nl.years <-22  # no. of no leap years in the data
start.date <- as.Date("1951/01/01")# start date in yyyy/mm/dd
end.date <- as.Date("2017/12/31") # end date in yyyy/mm/dd
date_seq <- seq(start.date,end.date,by = "days")

celltoind <-function(x){(x %% n.cols) + n.cols*(n.rows - (x %/% n.cols) - 1)}

*comment:Get the location coordinates
vilval<-cellFromXY(r,shp)

full_prec <- matrix(ncol = (366*l.years)+(365*nl.years),
                    nrow = length(shp), data = NA)
colval <- 0

for(num in seq(binfiles))

{
  myconn<-file(binfiles[num],"rb")
  #count=1
  while (length(a <- readBin(myconn, numeric(), n=n.cols*n.rows, size=4, endian = "little")) > 0 )
  {
    colval <- colval + 1
    a[a < (-90)]<- NA
    full_prec[,colval]<-a[celltoind(vilval)]
  }
  close(myconn)
}

full_prec<-as.data.frame(full_prec) # convert matrix to data frame

colnames(full_prec)<-date_seq # assign column names as date

prec_daily <- cbind(shp@data,full_prec) # this will assign shapefile dataframe

write_csv(prec_daily,"year_rf.csv") # write file to directory


texts:
I work with a raster set representing land values and I would like to automatically create polygons with random shape but with certain size and/or land value given a centroid point of each polygon. The centroid point will act as the initial raster cell (a pixel called as a seed cell) that will belong in that certain polygon and then a region growing algorithm will begin adding pixels around a centroid pixel until the constraint of size and/or land value of a polygon will be satisfied.
Is anybody who has any idea about an existing region growing algorithm? or at least which is the class of ArcObjects or Interfaces that I may work with, so as to read the value of each pixel, create a new polygon by adding new pixel to it?
I hope that the above make sense and somebody may help me.
Edit
Basically I intend to automatically create polygons with regular shapes and a series of constraints such as size, land value etc. using a genetic algorithm (GA). But first I need to create a set of random solutions in order to feed the GA.So, for instance on land block which is an area enclosed by roads I have already saying 4 centroid points. Each centroid point represents the approximate location of each new land parcel. In addition, each centroid is associated with attributes of each parcel i.e. size and land value. So, I want to begin by each centroid as a seed point of the region growing algorithm to starting creating a random shape for each parcel based on each centroid.
I hope that the above make sense. I look forward how can I create this region growing algorithm or if there is already one in VBA and ArcObjects.


texts: Point in Polygon using Rtree and Shapely

texts:
I am wondering what the correct syntax is for doing a point in polygon search using Shapely and RTree. I'm really interested in speeding the process up using the RTree index but am not having any luck getting it working with the "contains" function. See below.

First. I build my index.

idx = rtree.index.Index()
for feature in world:
    fid = int(feature['properties']['UN'])
    geom = shapely.geometry.asShape(feature['geometry'])
    idx.insert(fid, list(geom.bounds))

Then I check to see if the point is in the polygon using the contains function in Shapely/GEOS.

def point_inside_polygon_world(x, y):
    try:
        shapely.speedups.enable()
        for feature in world:
            #only work on features that match the country name
            shape = shapely.geometry.asShape(feature['geometry'])
            point = shapely.geometry.Point(float(x),float(y))
            check = shape.contains(point)

How would I check to see if the polygon from the index I created contains the point?

The contains function is not in the RTree docs so I was thinking it could potentially be something that's in there just not documented. Again, I am really looking to speed up the point in polygon lookup.


texts: Script that runs automatically from the QGIS Python Console when QGIS starts, using schedule tasks on Windows 10

texts:
I'm working with QGIS 2.8 and windows 10.
Is it possible that after open QGIS using the schedule tasks from windows 10, I could open the python console and run a script within QGIS programmatically.
How can I automate the following actions after open the QGIS:

open python console
open a script
run a script

Any help?

EDIT:

I've tried the two ways to open the python console, as suggested by Germán in this post "How to show Python console at QGIS program start" and it works. The problem occurs when I want to run my script. I tried to add my code these code lines as you suggest to define a startup.py file in my directory (C:\Users\USER\.qgis2\python), but as soon as I open QGIS I get this error:

"Traceback (most recent call last): File "", line 1, in File "C:/PROGRA~1/QGISWI~1/apps/qgis-ltr/./python\qgis\utils.py", line 478, in _import mod = _builtin_import(name, globals, locals, fromlist, level) File "C:/PROGRA~1/QGISWI~1/apps/qgis-ltr/./python\qgis\user.py", line 38, in execfile(startuppy, locals(), globals()) File "C:/Users/USER/.qgis2//python\startup.py", line 17, in layers =iface.legendInterface().layers() NameError: name 'iface' is not defined"


texts: Using RasterCalculator in ArcGIS Pro ModelBuilder gives ERROR 000539

texts:
As part of a larger model for iterating lake bathymetry interpolations, I'm essentially trying to establish the lake altitude as a new base level (Z=0) for the DEM surrounding the lake. (picture below is just one part of the full model)

My latest atempt at a workaround is trying to use RasterCalculator to "simply" subtract the aproximate lake altitude I've extracted from the local lake DEM (through Zonal statistics) using the equation pictured below. Where the "Value" variable is set as string.

However, when running it I get following error msg:
ERROR 000539: Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 5, in rcexec
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 4341, in Minus
    in_raster_or_constant2)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 4338, in Wrapper
    ["Minus", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError: ERROR 000732: Input Raster: Dataset 138.119995117 does not exist or is not supported
 Failed to execute (Raster Calculator).

I feel like there's probably a simple solution to this, but I really don't know where to go from here. From trying to look at the Functions code, I think it has to do with RasterCalculator not beeing able to use my "Value" variable. I've dabbled a bit with Python, but not to the extent to where I could solve this by just writing my own piece of code.


texts: What is topojson and how is it different than geojson?

texts:
This is a conceptual question.

What are the differences between these? and what are the pros and cons of these formats?


texts: Making definition query table not blue in ArcGIS Pro

texts:
I am new to ArcGIS Pro.

I figured out a definition query but it's making the attribute table blue and it hurts my eyes.   I would like definition query tables to be un-colored (white), but I want to retain colors for selecting from within tables using the normal selection process.

How do I do this?


texts: PyQGIS adding layer with auth config ID leaks into project

texts:
I add a layer using an authentification config ID using PyQGIS:

uri = QgsDataSourceUri()
uri.setConnection(HOST, PORT, DB, None, None, QgsDataSourceUri.SslPrefer, pg_auth_id)
QgsVectorLayer(uri.uri(), "mylayer", "postgres")

In the project, the credentials are written in clear text in the project file:

<datasource>dbname='pub' host=geodb.rootso.org port=5432 user='xxxxxx' password='xxxxxx'

Is it possible to avoid this? How?


texts: Adjust text font size in plotting autofitvariogram in R?

texts:
It should be quite simple to adjust the font size in plotting the outcome ofautofitVariogram using the automap package. The text fonts for the number of pairs are too big; the text fonts for the model, nugget... information are too small and not presentable. I was able to adjust the size of the points, but I failed to adjust the font of the label (or remove them) and the model information.

I've tried:

library(automap)
data(meuse)
coordinates(meuse) <- c("x","y")
data.auto.variogram <- autofitVariogram(formula = zinc ~1, input_data = meuse)
plot(data.auto.variogram, cex=1, par.settings=list(fontsize=list(text=5),par.main.text=list(cex=3), par.xlab.text=list(cex=2),par.ylab.text=list(cex=2)))

But that also reduces the size of the axis label and the information text.

We can use meuse as a reproducible example.

Can anyone give some solutions or recommendations please?


texts: ArcGIS Server 10 installation possible in Ubuntu 10.04?

texts: I have a machine with Ubuntu 10.04 LTS and want to install ArcGIS Server 10 in it. I find that the supported OS are RHEL and SUSE Linux in ESRI website and I also find in some postings that it is not possible in Ubuntu. Does anyone have any idea about it?


texts: Finding distance/length after exploding at vertices in QGIS

texts:
I have a shape-file consisting a road network. I wanted to explode the poly-lines at the vertices which I have manged to do. Prior to exploding I also used the function round($length,3) to find the distance of the poly-lines using the field calculator in the attribute table, and this gave the length of each poly-line.

I want to do the same after exploding to find the length of the exploded lines. But this time around I didn't get any values, instead got null values. Any idea how I can find the length of each exploded line?


texts: Different pixel values for a raster between QGIS and PostGIS

texts:
I have inserted a Landsat raster into PostGIS using the command :

raster2pgsql -I -C -s 28356 -M -t 100x100 /data/LS7_20150425_B40.tif public.red | psql -d postgres -U postgres

Using QGIS, I have opened both my original file and the raster from PostGIS using DB manager of QGIS.

According to QGIS the pixel values are not the same :

original file : values from 64 to 3687
raster from PostGIS : values from 10 to 8295

Note: if I export my raster from PostGIS to a new file using GDAL, the file is identical to the original one.

Any explanation about this difference ?


texts:
I am trying to run a few queries on spatial layers within a PostgreSQL/PostGIS database that I have connected to r=trough QGIS version 3.14. The queries themselves seem to run fine but I am unable t find a way of viewing the results in the QGIS map view. The 'load as new layer' option doesn't appear to work and reports 'an error creating the SQL layer'.
Can anyone advise on how to deal with this?


texts: Use if statement to judge in the Earth Engine

texts:
I hope to calculate such a formula through GEE：

But I don’t know how to write judgment statements and load this in GEE, Hope you can point out the error. Here is my code:
var getcsi = function (image) {
    if (nir.gt.ndwi) {
        var csi_1 = image.expression('sei - nir')
    }
    else {
        var csi_2 = image.expression('sei - ndwi')
    }
    return image;
};

https://code.earthengine.google.com/ca7f010a6ccd0b50bb32188bde5c6599


texts: Understanding tiles folder and file names

texts:
I have hard times finding the meaning of tiles names and folders. I don't think it's related to QGIS (it seems to be a standard thing) but when I export my map into a tile set the content is split into several folders.

The top folder is the zoom level. But inside I find a set of folders with numbers as names, and inside them a set of tiles with numbers as name. I understand the logic, folder is a "column" and every file is a new "row" of the column.

I imagine that these numbers define some coordinates in a way. But I can't find the logic behind it. For example here is the meta data exported by QGIS alongside the tiles :

"METADATA": {
 "map.coverage.center.longitude": 2.3287153244018484,
 "map.coverage.bottomRight.latitude": 48.83833969137705,
 "map.shortAttribution": "",
 "map.type": "overlay",
 "map.coverage.bottomRight.longitude": 2.343606948852528,
 "map.tileFormat": "png",
 "map.maxZoom": 18,
 "map.coverage.topLeft.latitude": 48.85906816414893,
 "map.coverage.center.latitude": 48.84870392776299,
 "map.shortName": "nc",
 "map.coverage.topLeft.longitude": 2.3138236999511683,
 "map.longAttribution": "",
 "map.minZoom": 16,
 "map.longDescription": "",
 "map.tileSideLength": 512
}

Is there a relation possible between the lat / long coordinates and folder / file names?


texts: How to use addImageQuery from leafem package with leafletProxy() in R shiny app

texts:
I'm writing a shiny app that generates a raster based on user input, and I would like to add this raster to a leaflet map and then have the raster values be displayed using addImageQuery from the leafem package.

So far I've only been able to use the addImageQuery function when I initially generate the map (see code block 1). However, when I first generate the map and then use leafletProxy to add the raster and image query, the raster is added to the map but the query does not work (see code block 2).

Code block 1 (addImageQuery works):

library(leaflet)
library(plainview)
library(leafem)
library(shiny)

ui <- fluidPage(
  fluidRow(
    leafletOutput("map")
  )
)

server <- function(input, output) {
  output$map <- renderLeaflet({
    leaflet() %>%
      addProviderTiles("OpenStreetMap") %>%
      addRasterImage(poppendorf[[1]], project = TRUE, group = "poppendorf",
                     layerId = "poppendorf") %>%
      addImageQuery(poppendorf[[1]], project = TRUE,
                            layerId = "poppendorf",position="bottomleft",prefix = "") %>%
      addLayersControl(overlayGroups = "poppendorf")
  })
}

shinyApp(ui, server)

Code block 2 (raster added on click, addImageQuery does not work):

library(leaflet)
library(plainview)
library(leafem)
library(shiny)

ui <- fluidPage(
  fluidRow(
    leafletOutput("map")
  )
)

server <- function(input, output) {
  output$map <- renderLeaflet({
    leaflet() %>%
      addProviderTiles("OpenStreetMap") %>%
      setView(11.2,49.73,12)
  })

  counter <- reactiveValues(
    counter = 0
  )

  observeEvent(input$map_click,{
    counter$counter <- counter$counter + 1
    if (counter$counter == 1) {
      leafletProxy("map") %>%
        addRasterImage(poppendorf[[1]], project = TRUE, group = "poppendorf",
                       layerId = "poppendorf") %>%
        addImageQuery(poppendorf[[1]], project = TRUE,
                      layerId = "poppendorf",position="bottomleft",prefix = "") %>%
        addLayersControl(overlayGroups = "poppendorf")

    }
  })
}

shinyApp(ui, server)

Does anyone know a way to get the addImageQuery function to work with a leafletProxy call?


texts: Geometry loss/mix up of a spatialite view joining two tables

texts:
I have two Spatialite tables:
- Sites (with geometry column) (ca. 5000 single sites)
- Samples (ca. 6800 samples)
To display the samples in QGIS I want to create a spatial view showing all samples at their respective locations from sites joined using the attribute field ID_site, however, the geometry is not joined (correctly) for all samples.

a) View Composer of spatialite_gui joins the tables using 'ID Site' and using rowid as the views_geometry_column - but always shows a "filtered view" displaying only the sites.
QGIS shows a "filtered view" restricted to the number of sites. All locations with more than one sample can not be accessed in QGIS via "Zoom to Feature/Flash Feature". Following the advice of Spatialite view appears 'filtered' in attribute table, won't load in canvas> Update SQL layer in QGIS with other unique(!) fields for all samples (e.g. id_sample or rowid from the sample table) do not change the result.

Variation of the join (FROM samples JOIN sites ON ... to FROM sites JOIN samples ON ...), variation the type of join (join, left join) as well as reassigning the id in QGIS via Right click layer Update SQL Layer do not result in the correct display of all samples.

b) by hand
Creating the View via SQLite and registering the view by hand (id_samples correctly registered in the view_geometry_columns table and updated via "Update SQL Layer" in QGIS as suggested here
Spatialite view appears 'filtered' in attribute table, won't load in canvas
 ) shows all 6800 samples in QGIS, however, their location is mixed up as the following images of a rectangular selection ares shows

CREATE VIEW view_join AS
    SELECT spls.ID_site,
           spls.ID_sample,
           sites.ID_site,
           sites.site_nom,
           sites.geometry
      FROM Samples AS spls
           JOIN
           Sites AS sites ON (spls.id_site=sites.id_site);
INSERT INTO views_geometry_columns (view_name, view_geometry, view_rowid,f_table_name,f_geometry_column,read_only)
VALUES ('view_join', 'geometry','id_sample','sites','geometry',0);

I am sure I am not the first one to join many samples to a few sites thus this complex problem for a standard routine took my by surprise - or I am missing something quite simple...
The problem is probably related to the joining or the adressing of the geometry column but having been through all variations I am at a loss and appreciate your help.


texts: PostGIS ST_AsPNG faster than ST_PixelAsCentroids

texts:
I'm hoping to export a PostGIS raster into a non-GIS enabled database. Ideally, I'd like to be able to have three columns in the output: latitude, longitude, raster_value. Easy enough.

Ideally, I'd accomplish this something like the following:

CREATE TABLE
    le_output
AS SELECT
    val,
    ST_AsText(geom) AS coord
FROM
    (
        SELECT
            (ST_PixelAsCentroids(rast, 1)).*
        FROM
            raster_table
    ) AS inside;

The problem is that this runs really slowly (hours). However, I can do something like this in a matter of 2 or 3 minutes:

SELECT
    ST_AsPNG(rast)
FROM
    raster_table;

However, if I go that route, now I have to do some post-processing in Python or some other language to open the PNG, process all the pixels, do the math against the raster grid and scale to find the lat/longs, etc. It's a huge pain I'd like to avoid.

Anybody have suggestions for a higher performance way to accomplish this? All I want is a list of raster pixel coordinates and their values, which seems like a trivial operation. I'm not sure why it's so darn fast with ST_AsPNG but so slow otherwise.

(For what it's worth, I've also tried joining the raster table to a list of every x/y coordinate location possible in the raster and using ST_Value. That's slow as molasses as well.)

Edit: As requested, the rasters are about 3500 x 3500 pixels, on average.


texts:
I am often creating hex grid layers in Leaflet which use circlemarker for styling.

I would like to shift / offset the shadow circle under each marker in order to give the layer some depth.

Which would be the best way to achieve this?

Could I possibly create a style function which creates two circlemarker for each feature with the lower circle being shifted and opaque.

I could just duplicate the layer. But how can I apply a shift to a layer or a circlemarker?

At the moment I just add a second layer with a lat/lon shift to the elements based on the zoom level.

var gridZoomShift = {
  12: 0.0004,
  13: 0.0002,
  14: 0.0001,
  15: 0.00005,
  16: 0.000025,
  17: 0.0000125
};

Not perfect but good enough for now.

I chose in the end the L.divIcon solution due to browser compatibility with filter:drop-shadow.

new L.geoJson(dots, {
  pointToLayer: function (feature, latlng) {
    return L.marker(latlng, divStyle(feature));
  }
}).addTo(map_geodata);

function divStyle(feature){
    var c = feature.properties.c;

    var s = c < circleRadiusArray[0] ? 4 :
        c < circleRadiusArray[1] ? 4.5 :
            c < circleRadiusArray[2] ? 5 :
                c < circleRadiusArray[3] ? 5.5 :
                    c < circleRadiusArray[4] ? 6 :
                        c < circleRadiusArray[5] ? 6.5 :
                            c < circleRadiusArray[6] ? 7 :
                                c < circleRadiusArray[7] ? 7.5 :
                                    c < circleRadiusArray[8] ? 8 :
                                        c < circleRadiusArray[9] ? 8.5 :
                                            9;

    return {
        icon: L.divIcon({
            className: 'circle',
            iconSize: [s, s]
        })
    };
}

The style for a bivariate map:

function divStyle(feature){
    var c = feature.properties.c;
    var v = feature.properties.v;

    var s = c < circleRadiusArray[0] ? 4 :
            c < circleRadiusArray[1] ? 4.5 :
            c < circleRadiusArray[2] ? 5 :
            c < circleRadiusArray[3] ? 5.5 :
            c < circleRadiusArray[4] ? 6 :
            c < circleRadiusArray[5] ? 6.5 :
            c < circleRadiusArray[6] ? 7 :
            c < circleRadiusArray[7] ? 7.5 :
                                       8.5;

    var circle_colour = v < circleColorArray[0] ? 'circle_c51b7d' :
                        v < circleColorArray[1] ? 'circle_de77ae' :
                        v < circleColorArray[2] ? 'circle_f1b6da' :
                        v < circleColorArray[3] ? 'circle_fde0ef' :
                        v < circleColorArray[4] ? 'circle_f7f7f7' :
                        v < circleColorArray[5] ? 'circle_e6f5d0' :
                        v < circleColorArray[6] ? 'circle_b8e186' :
                        v < circleColorArray[7] ? 'circle_b8e186' :
                                                  'circle_4d9221';

    return {
        icon: L.divIcon({
            className: circle_colour,
            iconSize: [s, s]
        })
    };
}


texts: How to create non-duplicate lines between all selected points in Spatialite?

texts:
There is a nice answer for the similar question for PostGIS.

I tried to adapt it for my case where I need to connect all points from a point layer, but unfortunately I get duplicate geometries: lines from point n to point m and from point m to point n. Here is the query I use. check_points is the point layer and cp_path is the layer where I want to store created lines.

WITH points1 AS (SELECT geometry FROM check_points),
     points2 AS (SELECT geometry FROM check_points)
INSERT INTO cp_path (geometry)
SELECT * FROM (SELECT DISTINCT MakeLine(points1.geometry, points2.geometry)
               FROM points1, points2
               WHERE points1.geometry != points2.geometry)

How can I avoid duplicated geometries?


texts: How can I force the blocksize when using gdal_rasterize?

texts:
I am rasterizing a shapefile of continents. I use the following command:

gdal_rasterize -burn 1  -ts 1440000 560000 -te -180 -60 180 80 -ot Byte -co COMPRESS=PACKBITS -co NUM_THREADS=ALL_CPUS -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 -co PREDICTOR=2 continent.shp continents.tif

You can see I specify the blocksize to be 256 x 256. Yet when I check the blocksize:

blah = gdal.Open(continents.tif)
blah.GetRasterBand(1).GetBlockSize()

Out[74]: [1440000, 256]

gdal_rasterize is overwriting my BLOCKXSIZE specification and instead using the number of columns. Is there any way to force the 256x256 blocksize?

I understand that larger x-blocksizes are more I/O efficient generally, but not for my specific application.


texts: Is it possible to access and manipulate the parts of the base topographic map from ArcGIS Pro?

texts:
I noticed the base topographic map that you start projects with has extensively detailed roads/streams/etc. for my area of interest and it's actually better than what I'm able to find elsewhere. This is the map in question. Does anyone know if there's a way to do one of the following:

Access/download a broken-up version of this map where I can access the individual vectors

Break up the topographic map directly to use said vectors

Find the specific sources they used to create this map

I'm expecting the answer to be "no it's proprietary" but I thought I'd check just in case.


texts: Optimal alignment of satellite images taken separately over the same geographic surface

texts:
I am given two satellite images of a location which has been subject to a natural disaster, one before the disaster and one after.  I am given bounding boxes for the images.The images have been aligned so that objects at ground level will appear at the same location in both images. For example, here is a picture of a shopping center about to be hit by Hurricane Harvey:

Here are the same bounding boxes drawn on the aligned post-disaster image:

It is visually evident that the bounding boxes are not well aligned with the post image, for objects which are not at ground level.  Per comment from Ian Turton below, this is due to different lookdown angle used for each image.

Suppose the pre and post images are both 1024x1024x3 RGB images. These are represented as NumPy arrays A and B. I would like to find a  transformation T(B) which best corrects for lookdown angle so that the bounding box is better aligned on B.  For example so that np.linalg.norm(A-T(B)) is minimized.   Is there a Python package that will do this for me, and how is the transformation expressed?  Assume that allowable transformation operations include rotation, translation, shear, and in general any kind of stretching or shrinking of the "fabric" of B so that it most closely lays on top of A.

This 2016 paper by Wang et al proposes a solution.  Prior to reading it, I imagined this kind of solution:

Project the 2D points into 3D space, for both pre and post
Infer the angle of attack and satellite height for both pre and post images
Do a perspective transformation from post angle/height to pre angle/height
Project the transformed 3D post image back to 2D

NOTE: The images above come from the xView 2 Challenge


texts: Create dynamic circle/polygon from specific lat/long using GeoTools

texts:
I have a dataset of lat/long that I would like to match on dynamically generated circle/polygon from specific point with different radius.

I have already completed how to match lat/long with custom defined boundaries in a shapefile using GeoTools.

Now, instead of using boundaries defined in a custom shapefile, I would like to do following steps

Find centroid of the boundary (polygon)
Using centroid lat/long, create multiple circles/polygon with radius 5, 10 and 15 miles like  5 miles = Polygon1,   10 miles = Polygon2  15 miles = Polygon3
Match dataset of lat/long to see how many points lie in
    Polygon1, Polygon2 and Polygon3.


texts: How can I open a jpg format with QGIS 2.8?

texts: I'm having some difficulties opening the jpg format at QGIS 2.8. It opens, but nothing appears.


texts: How to find the distance between start of line and start of raster with arcpy/python?

texts:
I'am creating a python script tool for ArcGIS Desktop 10.1 to extract profiles. I’m using stack profile to extract values to a table but stack profile only outputs values where the line and raster overlaps.

My solution to create profiles that will start at the beginning of my line is to find the distance between the start of the line and the start of the raster and then offset the x-axis. But how do I find this distance (marked with blue in picture below) with arcpy/python?

Update: the solution needs to work with floating point raster dataset


texts: Convert MODIS_Sinusoidal into WGS84 in Google Earth Engine

texts:
I'm trying to transfer the Modis_Sinusoidal projection to the WGS84 projection in the Google Earth Engine. But I failed. I want to know how to solve it.
var bp = ee.FeatureCollection("users/miaohu828/bankment_polygon"),
    modisLandcover = ee.ImageCollection("MODIS/006/MCD12Q1");
Map.setCenter(116.05380098858659, 30.63744097417238, 6);
var filtered2018 = modisLandcover.filter(
ee.Filter.date('2018-01-01', '2018-12-31'))
var landcover2018 = ee.Image(filtered2018.first())
var classified2018 = landcover2018.select('LC_Type1')
// The input image has a SR-ORG:6974 (sinusoidal) projection.
// Operations *before* the reproject call will be done in the projection
// specified by reproject().  The output results in another reprojection.
var reprojected = classified2018
    .unitScale(-2000, 10000)
    .reproject('EPSG:4326', null, 500);
    var visParams = {min: 1, max: 17};


texts: What tools do you recommend for this kind of web map?

texts:
I have one tourist map (ski resort) with buildings and roads in .jpg format. It's problem because north is down and south is up (it is that kind of drawn map with drawn mountain up on map where is actually south).

I need to put it on the website and also put some pins on it.

What tools and ways do you recommend me? It needs to be done fast so please give me some simple solutions if possible.


texts:
I have a script in Python that includes a series of Spatial Analyst tools.  I need to reclassify a raster and have tried using Remap which fails and have also tried a Con statement.

The Con statement is:

outCon = Con(IsNull("median141516_nn_test"),0,Con("median141516_nn_test" <= 5.0,(("median141516_nn_test"-0.0)/(5.0-0.0)*(1.25-0.0)),Con("median141516_nn_test" <= 10.0,((("median141516_nn_test"-5.0)/(10.0-5.0)*(2.5-1.25))+1.25),Con("median141516_nn_test" <= 15.0,((("median141516_nn_test"-10.0)/(15.0-10.0)*(3.75-2.5))+2.5),Con("median141516_nn_test" <= 20.0,((("median141516_nn_test"-15.0)/(20.0-15.0)*(5.0-3.75))+3.75),Con("median141516_nn_test" <= 25.0,((("median141516_nn_test"-20.0)/(25.0-20.0)*(6.25-5.0))+5.0),Con("median141516_nn_test" <= 30.0,((("median141516_nn_test"-25.0)/(30.0-25.0)*(7.5-6.25))+6.25),Con("median141516_nn_test" <= 65.0,((("median141516_nn_test"-30.0)/(65.0-30.0)*(8.75-7.5))+7.5),Con("median141516_nn_test" > 65.0,((("median141516_nn_test"-65.0)/(100.0-65.0)*(10.0-8.75))+8.75))))))))))

# Save the output
outCon.save("rcl_141516_test")

And the error is:

Runtime error
Traceback (most recent call last):
  File "", line 1, in
TypeError: unsupported operand type(s) for -: 'str' and 'float'

The raster is a floating point.  The only thing I can think is that Python thinks the raster is a string, but I have saved a raster after each tool (similar to outCon.save) and have not received this error

I am new to writing multi tool scripts.


texts: I am using online Google Maps links for imagery in QGIS. If their is any network error I am unable to work. So please suggest me if their is any offline maps for QGIS?


texts:
I need to interpolate the 3" srtm data for a game I'm developing (it's like Minecraft - voxel style). I downloaded the srtm data from viewfinderpanoramas.org. My aim is to get the height for every meter, so I'm tiling the world into 16x16m chunks so I have to perform 256 interpolations per chunk. For every block in a chunk I'm doing the interpolations like that:

Calculate lat/lon for the block my chunk with the help of this post
Calculate lat/lon for the 4 points q11 q12 q21 q22 which I need for my bilinear interpolation. You can see this calculation on the following picture. I'm add/sub 45m to each block's lat/lon because I want to find the nearest grid-point for each q11 q12 q21 q22 in my srtm-file. (the picture shows the calculation for the block in the top-left corner)

Then I'm searching the height data for q11 q12 q21 q22 in my srtm-file with the help of this post
I'm performing my bilinear interpolation

The problem is, that my interpolated data looks like a heap of ruins but you can recognize the approximate terrain.

Can you tell me some hints or mistakes concerning my method?

Especially hints how I can get the 4 nearest grid points for the interpolation.


texts: Zoom in, zoom out, and pan tools cannot be toggled off in QGIS?

texts:
When I select the Pan Map, zoom in, and zoom out tools in the QGIS toolbar (version 2.18.16), I cannot toggle it off after using it. In efforts to toggle a tool from the toolbar off after using it, I have tried clicking the tool again and pressing escape but none works.

How do I toggle off a tool like pan or zoom in QGIS?


texts: Accessing "many to many" relationship class using ArcGIS Identify?

texts:
I am trying to create a many to many relationship class in ArcGIS 10.0. When Identifying the features though, I don´t see the associated table values.

Here is the Class schema:

Destination Table name = "TypeDef" ( a table containing a list of Type definitions)

TypeID* | Typ_Name

Origin Table name = "Points" ( A point geometry)

ElementID* | Element_Name

Association Table name ="ElementTypes" (join table with ElementIDs and TypeIDs)

ElementID* | TypeID*

My ManyToMany Relationship class is defined as:

SourceTable
PK: ElementID
FK: ElementID

DestinationTable
PK: TypID
FK: TypeID

However, when I use the Identify Tool on a point, only the Relationship name is shown and not the attributes. What am I doing wrong?


texts: Create an optimised path between nodes avoiding polygons entities

texts:
I have a point shapefile representing seaports. I have the linear cost of shipping for every seaport pairing (specific shipping cost to go from one point to another point). I want to find a way of creating shipping routes (lines) between all nodes.

Then I will be able to give, for every line, a specific linear cost ($/km), so I can calculate, if I want to go from seaport A to Seaport Z, which path is cost optimised through pgRouting.

On top of that, as a constraint, I want the shipping to avoid land (I have a shapefile of land polygons). That's the real difficulty, because that's way more difficult than a "shortest path" method and I can't figure out any method to deal with this difficulty as a novice.
I hope I  was clear enough.


texts: Are tab file in QGIS the same as file in MapInfo?

texts: So I received a bundle of files under the same name with the file formats of dat.;map.;id.;and tab. After searching on the web, i am assumming that this type of file is most probable to be produced in MapInfo. However, I only have a QGIS and was only able to open tab. format. I was wondering if the information in the tab. file will the be same that was produced from MapInfo? or is it a file produced from a complete different GIS system.


texts: QGIS 3.18.1 Print layout grids are incorrect, cannot show them in DMS format

texts:
Normally, in ArcMap Print layout, enabling the grids and specifying the intervals would show the curved grids as in the following map:

But in QGIS, enabling the grids (1) does not have an option of showing the coordinates in DMS, (2) reveals wrong coordinates, (3) is not curved like in the previous map:

I assumed it was an error with the projections and made sure that the map was assigned the same CRS (WGS84 - UTM Zone 48N) as the project, but the last photo was still the result.
I saw a "solution" that changing the projection to EPSG:4326 WGS84 would supposedly allow for DMS units to be shown, but changing the projection to WGS84 returned an ugly output:

where the coordinates were still incorrect.
How would one go about (1) displaying the correct grid labels as DMS, (2) show the gridlines in a curved fashion on a map projected to UTM on QGIS 3.18.1?
Grid settings:


texts: Create a esri layer file through web map

texts: Is there some way to create a layer file through web mapping application.  If I expose data through a map service is there a way to create a .lyr file?


texts: Distance to features in R raster with a buffer?

texts:
Is there a method to calculate the distance to features in a raster, but only within some buffer distance? At the moment I'm just using distance(rast), and then setting all the values greater than the buffer  to be equal to the buffer: rast[rast > buffer] = buffer
Here's an example:
 r <- raster(xmn=1035792, xmx= 1116792, ymn=825303.6, ymx=937803.6, resolution = 12.5,crs = "+proj=aea +lat_1=42.122774 +lat_2=49.01518 +lat_0=45.568977 +lon_0=-84.455955 +x_0=1000000 +y_0=1000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")
r <- setValues(r, 0)
x <- c(1088419, 1093067, 1090190, 1087977, 1088419)
y <- c(895030.8, 894366.9, 892153.9, 893703.0, 895030.8)
poly <- SpatialPolygons(list(Polygons(list(Polygon(data.frame(x,y))), ID = 1)))
r <- mask(r, poly)
rD <- distance(r)
rD[rD > buff] = buff

This is just a small, very simple example, but some of the rasters I'm using are extremely large, and can be much more complicated than just this little polygon (for example road networks over whole countries), and the distance function can be extremely slow when this is the case. So, I'm wondering if there's a better method, that will only calculate the distance up to a maximum buffer and hopefully save on computation time? I've spent a long time searching and haven't found anything useful, so hopefully someone can help me!
Thanks in advance.


texts: My layer won't appear in QGIS using PostgreSQL 9.1

texts:
I am working with a "cities" layer (polygons) in a PostgreSQL 9.1 database.  When I use the following query in PostGIS:

select st_buffer(st_union(st_centroid(cities.geom)),1500) geom from "my_layer" cities

I have no problem with displaying the layer in QGIS.

When I  use this syntax:

select st_union(st_buffer(st_centroid(cities.geom),1500)) geom from "my_layer" cities

No error is generated in PostGIS, the generated table has one element (as expected) the table can be opened in qgis, its extent "seems" OK but nothing appears on the screen ...

I checked objects with a ST_IsValid query and found no error.  I tried adding an id but nothing changed.  I tried these two queries in PostgreSQL 9.5 and had no problem at all.

What could I be doing wrong, why do I have this behavior using PostgreSQL  9.1, and what could be the workaround?

PS : PostGIS version is 1.5.3 in the first case (PostgreSQL 9.1) and 2.2.1 in the second case (I suppose it must be important ...)


texts: "Computation timed out" when trying to get the size of an ImageCollection in Google Earth Engine python installation

texts:
I'm trying to get the size of an ImageCollection but get "computation timed out" while processing the getInfo() request. I'm using the python jupyter installation and I'm unsure how to make a callback as suggested in the javascript tutorial.

The request is as follows:

area_of_interest_sen = ee.ImageCollection('COPERNICUS/S2').filterDate('2017-06-01', '2018-06-01').filterBounds(jsonf)\
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
s = area_of_interest_sen.size()
print(s.getInfo())

Where jsonf is a geometry from a file.

The error goes like this:

EEExceptionTraceback (most recent call last)
<ipython-input-215-e097e3c89cef> in <module>()
     58 s = area_of_interest_sen.size()
     59
---> 60 print(s.getInfo())
     61
     62 area_of_interest_sen = area_of_interest_sen.median().select(bands)

/src/earthengine-api/python/ee/computedobject.py in getInfo(self)
     93       The object can evaluate to anything.
     94     """
---> 95     return data.getValue({'json': self.serialize()})
     96
     97   def encode(self, encoder):

/src/earthengine-api/python/ee/data.py in getValue(params)
    253   """
    254   params['json_format'] = 'v2'
--> 255   return send_('/value', params)
    256
    257

/src/earthengine-api/python/ee/data.py in send_(path, params, opt_method, opt_raw)
    795       raise ee_exception.EEException('Invalid JSON: %s' % content)
    796     if 'error' in json_content:
--> 797       raise ee_exception.EEException(json_content['error']['message'])
    798     if 'data' not in content:
    799       raise ee_exception.EEException('Malformed response: ' + str(content))

EEException: Computation timed out.

PS: other than this, the query works. I get the images and I'm able to process through them.


texts: Predefined CRS and Custom CRS problem with Borneo projection

texts:
I loaded a shapefile layer into QGIS and assigned to CRS of EPSG:29873 (Timbalai 1948/ RSO Borneo (m))

I copied the Proj4 parameters from EPSG:29873 (above) and create a custom user defined coordinate systems (as below), later I intend to add in false easting and northing of 2 mil / 5 mil.

+proj=omerc +lat_0=4 +lonc=115 +alpha=53.31582047222222 +k=0.99984 +x_0=590476.87 +y_0=442857.65 +ellps=evrstSS +towgs84=-533.4,669.2,-52.5,0,0,4.28,9.4 +units=m +no_defs

Before modification (adding false easting and northing), in QGIS, I inserted a shp layer (same source) twice but assigned them to two different CRS for testing (one to EPSG:29873, and the other one to user defined USER:100001, Basemap using Bing Aerial with project CSR to EPSG:3857 WGS84/Pseudo Mercator). the result displayed quite surprisingly. the EPSG:298873 CRS locate the shp file correctly but NOT the user defined custom CRS (300+++ Km toward north from the correct one)! In fact they should displayed at the same location right, seeing that they are using the same Proj4 parameters!!!

Did I create the custom CRS wrongly?


texts: Create a common scale for different raster layers in GRASS, SAGA or QGIS

texts:
I have a number of layers that represent the temperature variability of the water in a lake.

They have been generated  trough interpolation and have different values as they represent temperatures from different seasons.

Everything looks fine in the maps except that they have similar colors although they represent different values. For example on the March-May map I have red for 10 degrees C, while for June-September the same red color is used for 20 C.

Ideally I should be able to extend the temperature (value) range for all of them to be the same. So if I look at the spring map I should seem more blue (suggesting colder water)  while in the summer map more red regions.

I played with the d.legend and I can modify the range but only as a subset i.e if for spring I have 0-12 C I can reduce it to 0-8 C, but I cannot extend it to 0-20.

So now the question: How can I unify the legend ranges for all these layers so that a specific color represents the same value on all the maps?

I looked at r.recode and r.rescale, but they seem to take the same color range and change only the number corresponding to that color. I want to change the color for the value and not the opposite. In addition the rescaled scalebar is discontinued instead of a nice gradient.


texts: I have shapefile for my country and I want to view this map in my web site (browser) in my LAN network how to do that, is it possible?


texts:
I have a several hundred points that I am showing on a map. At the moment, my method of adding the points seems inefficient. I create an array using the following approach:

var features = [new OpenLayers.Feature.Vector(OpenLayers.Geometry.fromWKT("POINT(-4.2522971 45.3447353)").transform(new OpenLayers.Projection("EPSG:4326"), map.getProjectionObject()),{"Name":"aaa"}),
new OpenLayers.Feature.Vector(OpenLayers.Geometry.fromWKT("POINT(-4.259215 45.344827)").transform(new OpenLayers.Projection("EPSG:4326"), map.getProjectionObject()),{"Name":"bbb"})]

and then add this array to a map-layer:

var vectors = new OpenLayers.Layer.Vector("vector", {isBaseLayer: false, projection:"EPSG:4326", styleMap:myStyle});

vectors.addFeatures(features);
map.addLayer(vectors);

Is there are a more efficient way of displaying the points? I also want to be able to dynamically change the radius of each point separately. I'm not sure if that is something that I should consider when constructing the point. For example, I can choose the radius when constructing, but I do not know how to update it:

new OpenLayers.Feature.Vector(OpenLayers.Geometry.fromWKT("POINT(-4.2522971 45.3447353)").transform(new OpenLayers.Projection("EPSG:4326"), map.getProjectionObject()),{"Name":"aaa"}, {radius: 10})


texts:
I am trying to get the polygon data from the following leaflet web map http://www.cityfibre.com/rollout/.
I am not really sure what kind of format I should be looking for (geojson, xml, etc..) but I tried the following with no luck.
I go to inspect the element -> network tab -> XHR , I then refresh the page and check the response but can't find anything related.
Any ideas?


texts: Loading collada into 3D view of QGIS using Qgis2ThreeJs Plugin?

texts:
I have been looking at documentation and other posts but cannot find a way to load a 3D model into QGIS.

I Installed the qgis2ThreeJs Plugin.

I have tried to load a collada file in a temporary scratch layer. The Geometry type was Point.

After adding the layer I launch the threejs plugin.

Then I can go into the point layer properties and set the Shape to "3D Model". Then you can select a model file. After closing the properties, I do not see any model in the 3D view.


texts:
I have the following issue when using the Create Database View tool in ArcMap.
My query runs fine in Microsoft SSMS and Visual Code (with mssql extension) against the enterprise geodatabase.
But for some reason, when I'm trying to create the view using Create Database View in ArcMap, I get the following error.
Traceback (most recent call last):

  File "D:\Beni\OT\Databases\Actualizare\Tools\Scripts\10_Creare_view.py", line 39, in <module>

    arcpy.CreateDatabaseView_management(in_db,numeView,defView)

  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\management.py", line 4318, in CreateDatabaseView

    raise e

ExecuteError: ERROR 999999: Error executing function.

[[Microsoft][ODBC Driver 11 for SQL Server][SQL Server]Incorrect syntax near '('.]

DBMS table not found

Failed to execute (CreateDatabaseView).

ERROR 999999: Error executing function.

[[Microsoft][ODBC Driver 11 for SQL Server][SQL Server]Incorrect syntax near '('.]

DBMS table not found

Failed to execute (CreateDatabaseView).

I can successfully connect to the Enterprise Geodatabase in ArcMap and create other views using the same tool, on the same geodatabase.

I am able to create the view via Catalog window in ArcMap, using right-click on geodatabase connection - New - View. It works flawlessly this way and it doesn't complain about syntax.

I could create the views manually, via the above-mentioned method, but I am trying to automate the operation via the Create Database View tool. This is not the only query that fails with this tool (I have about 400 of them which work and about 70 which fail).

I checked the syntax a thousand times and I still can't figure it out. Maybe someone here could help me understand the issue. My only conclusion is that this could be a bug in the Create Database View tool, with certain types of queries.
The SQL query looks like this:
SELECT
    CAST (ROW_NUMBER() OVER (ORDER BY ADM.DENUMIRE) AS INT) AS OBJECTID,
    ADM.DENUMIRE,
    CAST(ADM.LOC_COMPONENTE_ALE_MUN_ORAS AS DECIMAL) / CAST(iif(AGR.ANUL_2014 = 0, -1, AGR.ANUL_2014) AS DECIMAL) AS DENSITATEA_ASEZARILOR_URBANE,
    ADM.LOC_COMPONENTE_ALE_MUN_ORAS,
    AGR.ANUL_2014 AS SUPRAFATA,
    ADM.SHAPE
FROM
(SELECT
    JUD.OBJECTID,
    JUD.DEN_JUD AS DENUMIRE,
    ADM.LOC_COMPONENTE_ALE_MUN_ORAS,
    JUD.SHAPE
FROM
(SELECT
    JUDET,
    SUM(LOC_COMPONENTE_ALE_MUN_ORAS) AS LOC_COMPONENTE_ALE_MUN_ORAS
FROM gisdata_owner.ADM010 AS ADM
GROUP BY JUDET) AS ADM
RIGHT OUTER JOIN
gisdata_owner.JUDETESIMP AS JUD
ON JUD.DEN_JUD = ADM.JUDET) AS ADM
INNER JOIN
(SELECT
    AGR.MACROREGIUNI_REGIUNI_DE_DEZVO AS DENUMIRE,
    AGR.ANUL_2014
FROM gisdata_owner.AGR101A AS AGR
WHERE AGR.FORME_DE_PROPRIETATE = 'Total'
AND AGR.MODUL_DE_FOLOSINTA_A_FONDULUI = 'Total') AS AGR
ON AGR.DENUMIRE = ADM.DENUMIRE

Setup info
Windows 10 Pro v. 2004 (build 19041.572)
ArcGIS Desktop Advanced / ArcMap 10.5.0.6491
SQL Server 2014 Developer Edition (12.0.6118.4)


texts: What's the best way to display thousands of polygons in Openlayers?

texts:
I am trying to create an Openlayers map that will display thousands of polygons. And when the user clicks on each, I want some of the metadata to display. Working with this large of a dataset seems like it will be challenging. But I am trying to figure out if there is a standard way to go about this. At this point, I am think I should either

display my polygons with WMS and use getfeatureinfo
Display with WFS, Load with Strategy.BBOX, and use GetFeature to get meatadata.

Is there a logical way I am missing to load all of this data into a map and have a good speed. I was looking at this page and trying to figure out what they did: http://protectedplanet.net/


texts:
I use QGIS 3.6 on Windows 10 and I want to reset it.
I tried this through Setting-Options-System-Reset user Interface to default settings. I uninstalled and installed again, but still there are some saved features.


texts: Changing Styling on Individual Leaflet Layer Control Items

texts:
How would I go about changing the styling for Individual entries in Leaflets layer control?

For example I would like to change the colour of "Layer 1" to red and the colour of "Layer 2" to blue.

Would this be done through the Javascript or the CSS?


texts: How can we generate GeoJSON from SQL data with latlong

texts: Is there any way we can create GeoJSON from SQL server table with Latitude & Longitude Data?


texts: How to show some labels only at certain scales?

texts:
I am using QGIS 2.14 and I want to set my contour line lables so that my minor labels (every half a metre) are only shown at a scale under 1:1000, while having all my major (every metre) labels always showing. I assume you can do this in the expression builder - but I am not sure where to find the 'Maximum' function name.

So far I have:

Case
When "Elevation" Like '%.5' THEN ????
When "Elevation Like '%.0' THEN ???
END


texts: Add WMS from GeoServer to Openlayers with OSM baselayer

texts:
I am new to JavaScript/OpenLayers.  I would like to add a WMS layer (BLD_504) from dwasa (workspace) locally. The code is as follows:

<!DOCTYPE HTML>

<title> Dhaka WASA </title>

When I zoom in I see the following:

The data looks good in GeoServer:


texts:
Update: The bug has been fixed in the ArcGIS 10.4 release

I am using ArcGIS 10.2.2 to determine zonal statistics for a number of zones. If there is any NoData in the value raster, I want the zone results to be "NoData", precisely as advertised by the tools description. This tool description states:

  DATA — Within any particular zone, only cells that have a value in the input Value raster will be used in determining the output value for that zone. NoData cells in the Value raster will be ignored in the statistic calculation.

  NODATA — Within any particular zone, if any NoData cells exist in the
  Value raster, it is deemed that there is insufficient information to
  perform statistical calculations for all the cells in that zone;
  therefore, the entire zone will receive the NoData value on the output
  raster.

Please have a look at my setup in this picture:

I am using the NODATA option with a value raster that has one NoData pixel, and therefore expect the resulting zone value (zone 61154) to be 'NoData'. Instead, I get a value of 12.74 (rounded to 13 in the image), which confuses me on two levels: First, I expected 'NoData', and second, the resulting value of 12.74 is mathematically impossible, because the mean cannot be larger than the maximum value in the value raster, which is 10 in this case.

If I am using the DATA option, I get a value of about 9.1, which makes sense. We tested this on different datasets, computers, and ArcGIS versions.

What am I missing here?

Edit / Additional comment: I just noticed that the 'Count' attribute is also wrong for that particular zone. There are indeed 421 cells in that zone, but the tool only counted 297. Calculating 421 minus 297 results in 124 - oddly enough, this is the "position" where the NoData pixel is located, if one counts the pixels from upper left to lower right in the zone. The tool might be getting the cell count wrong (too low), which might explain the increase of the average.

Edit: Here is a link to the data I am using.

Edit: Dan Patterson and I did some further debugging here at the ESRI forum.


texts: Incorrect time displaying in ArcGIS API for JavaScript popup template?

texts:
I'm creating a web map that is updated a few times a day showing voltage readings from meters on our distribution network.

This layer has a field called ReadTime, and the data type is Date.

The date is formatted as: 0/00/0000 00:00:00 AM

So, I'm using the shortDateLongTime format as shown here.

However, when I click on a meter the date comes in fine, but the time is wrong. For example, a value for one of the meters in the attribute table from ArcMap is 6/17/2015 8:00:14 AM, but in the web map it's displaying as 6/17/2015, 1:00:14 AM.

I'm assuming I'm missing a setting somewhere, but I haven't been able to find out where.

var popupTemplate = new PopupTemplate({
        title: "{VoltageReadings_LOCATION}",
        fieldInfos: [{
            fieldName: "VoltageReadings_READ_TIME",
            visible: true,
            label: "Read Time",
            format: {
                dateFormat: 'shortDateLongTime'
            }
        }],
        showAttachments: true
    });

Update to show what worked:

var popupTemplate = new PopupTemplate({
        title: "{VoltageReadings_LOCATION}",
        fieldInfos: [{
            fieldName: "VoltageReadings_READ_TIME",
            visible: true,
            label: "Read Time",
            format: {
                dateFormat: 'shortDateLongTime'
            }
        }],
        showAttachments: true
    }, {
        utcOffset: 0
    });


texts:
I am trying to convert a PostgreSQL table to CSV using ogr2ogr. Why is my command failing?

ogr2ogr -f CSV output.csv PG:"host=myhost port=5432 user='postgres' password='mypassword' dbname='mydb'" "myschema.mytable"

This is the error I get:

ERROR 1: PQconnectdb failed.
invalid connection option ""host"

FAILURE:
Unable to open datasource `PG:"host=myhost' with the following drivers.

All of the drivers are listed, including PostgreSQL.

I have read that syntax and placement of quotations are important, but nothing I try seems to work.

I am trying to run this command in a Ubuntu Docker container. I installed gdal like this:

apt-get install gdal-bin

Version info:

# ogr2ogr --version
  GDAL 1.11.3, released 2015/09/16

I have verified that I can connect to the database using other similar method. For example, when I run ogrinfo, I connect and can see a full list of spatial tables, including the one from which I am trying to export.


texts: Layers getting delinked every time project closes

texts:
I have a QGIS file with about 10 layers (shape files, web and a joined csv). Every time I close the file and open it, all layers are unlinked and the 'Handle Bad Layers' dialog appears. Needless to say, its tedious to relink every time and this behavior was not there earlier.

I've not changed the location of any of the layers. In fact, I tried moving them all to one location and then relinked them, but after I shut and opened the file again, the behavior did not change and the bad layers dialog came up again.

What could be going wrong? How could I correct this?

QGIS 2.18.16 on Ubuntu 16.04


texts:
I am creating a plugin for Qgis 2.6 and when my plugin loads up i want to disable the visibility of some of the tool bar. I was able to set the visibility for various tool bars like

    self.iface.digitizeToolBar().setVisible(True)
    self.iface.mapNavToolToolBar().setVisible(True)
    self.iface.advancedDigitizeToolBar().setVisible(True)
    self.iface.pluginToolBar().setVisible(True)
    self.iface.layerToolBar().setVisible(True)

    self.iface.databaseToolBar().setVisible(False)
    self.iface.fileToolBar().setVisible(False)
    self.iface.helpToolBar().setVisible(False)
    self.iface.vectorToolBar().setVisible(False)
    self.iface.rasterToolBar().setVisible(False)
    self.iface.webToolBar().setVisible(False)
    self.iface.attributesToolBar().setVisible(False)

Apart from this tool bars I  have other tool bars like Manage Layer, Grass, Label but I am not able to get the reference for this tool bar from the iface to set the visibility.

http://doc.qgis.org/api/classQgisInterface.html#a882199db87f1b58bd5fc9b3d5db00d28
This documentation does not expose this tool bars.


texts: Where can I find a python GDAL 1.8.0 package?

texts: I have python 2.7.4 installed on my RHEL 5 machine along with gdal-devel-1.8.0 and gdal-1.8.0. I tried to install the GDAL 1.9.1 package available on python's page (https://pypi.python.org/pypi/GDAL/1.9.1) but got errors that are related to trying to install v1.9.1 when my gdal libraries are 1.8.0 (see Error installing GDAL 1.9.0 python package). I can't find any other versions of python GDAL bindings. Can anyone help me please?


texts: Leica Zeno Office on ArcGIS 10.2.2 Easy In Issue

texts:
I have had Leica Zeno Office 3.2 installed on ArcGIS 10.2.2 for the last 3 months with no problems but for some reason I can't bring the field data back in with easy in.

See images below for my process:


texts: QGIS Plugins & Tools that replicate Select ArcPy Tools

texts:
I developed a complex python application for road generation from GPS Devices that uses a few Arcpy Tools.

I want to try and replicate this Applications functionality using a freeware like PyQGIS instead of Arcpy, but in order to do this I need to find tools & plugins from QGIS that mimic the following tools:

arcpy.da.SearchCursor
arcpy.TableToTable_conversion
arcpy.JoinField_management

arcpy.MakeQueryTable_management
arcpy.Delete_management
arcpy.AddField_management
arcpy.CalculateField_management
arcpy.MakeXYEventLayer_management
arcpy.CopyFeatures_management
arcpy.PointsToLine_management,  QGIS Points to Path: http://spatialgalaxy.net/2011/12/30/qgis-plugin-of-the-week-points-to-paths/
arcpy.Intersect_analysis
arcpy.SpatialJoin_analysis
arcpy.MakeFeatureLayer_management
arcpy.SelectLayerByAttribute_management
arcpy.LayerToKML_conversion
arcpy.ApplySymbologyFromLayer_management

I am also researching this and will post the matching tools in the question as I find them to help anyone who is having similar issues, once all Arcpy tools have a matching PyQGIS tool or plugin I will accept an answer with them all in it, If an answer has at least 1 matching tool I will upvote, but no answer will be accepted until each arcpy tool has a matching PyQGIS


texts: Reclassify in QGIS using Raster Calculator?

texts: I have downloaded data from the USGS to create a Map for a malaria risk in east africa but my question is only about the actuall Input in the Raster Calculator. I have a RFE Tiff that shows the Rainfall in the Area, the highest rainfall is within the value of 76 to 150 so I used the following expression in the Raster Calculator to get an Output with only these values: "ea06063_ea3606mod@1" >= 76 AND "ea06063_ea3606mod@1" <= 150. Unfortunately all I got as an Output was a black screen... Anybody knows what could be the mistake??


texts: Results in ArcGIS search window do not bring up tool when clicked

texts:
I have been unable to use the links resulting from a tool search to open a tool. When I click the link for the tool nothing happens - no tool is retrieved.

This problem is only occurring with tools and not maps, data, and images. Also, this problem occurs with every tool I have tried to open from the search window. Tools open regularly from toolbox window.

What do I need to do to restore functionality to the search window?


texts: Remove terrainProvider from scene (ol3-cesium)

texts:
I'm using ol3-cesium for creating 3D globe and cesium terrain provider to show terrain data on globe its working fine.

this.scene.terrainProvider = new Cesium.CesiumTerrainProvider({

url : "http://" + window.location.host + "/tilesets/world",

requestVertexNormals : true

});

But on click of a button i want to hide the terrain data or remove terrain provider from scene object, which ever is possible.


texts: How to join values from a R dataset to a shapefile of polygons I've created in QGIS?

texts:
I have a dataset where I have count data of sea turtles associated to different areas, for example in area 4 the divers saw 20 turtles.
I want to display the column of counts on a shapefile I've created in QGIS to make something like a heat map or a graduated scale color.
This shapefile consists on polygons I've designed and in the table of the attributes, the ID's of the polygons are the number of each area.
I've imported this shapefile into R and I've been trying to use geom_map but I don't know what to put on the aes. It asks me for x, y and id


texts:
I'm trying to iterate through a GeoJSON file (below) and eventually populate an array with a few attributes within "properties". I keep coming up short working through various examples I can find. Do I need to embed another $.each() to get to properties. If it isn't obvious I'm new to this and I've hit the wall.

What I have so far:

$(document).ready(function () {
            $.getJSON('testpoint.geojson', function (data) {
                var items = [];
                $.each(data.features, function (key, val) {
                    items.push('
' + val + '
');
                });
                $('
', {
                    'class':'my-new-list',
                    html:items.join('')
                }).appendTo('body');
            });
        });

And my json is as follows:

{
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "properties": {
                "gid": 1,
                "PID": "9-18-3",
                "BCT": "BCT126",
                "OWNER_TYPE": "A",
                "LOCNO": 0,
                "LOCEXT": "",
                "STREET": "CROSBY LANE",
                "ACQUIRED": "5/7/2010",
                "GRANTOR": "John A. SPARGO",
                "UPLAND": 0,
                "WETLAND": 3.96,
                "TOTAL": 3.96,
                "HABITAT": "salt marsh"
            },
            "geometry": {
                "type": "Point",
                "coordinates": [
                    -70.03209,
                    41.78278
                ]
            }
        }
    ]
}


texts: Missing qgis icon resources for symbology dialog

texts:
I'm trying to call QgsRendererV2PropertiesDialog from a standalone Qgis python script. It mostly works, except there's some "unknown error" messages related to SVG files when I trigger the dialog:

Cannot open file ':/images/themes/default/sort.svg', because: Unknown error
Cannot open file ':/images/themes/default/sort.svg', because: Unknown error
Cannot open file ':/images/themes/default/symbologyAdd.svg', because: Unknown error
Cannot open file ':/images/themes/default/symbologyAdd.svg', because: Unknown error
[etc]

and the relevant buttons that should have those icons are blank. The dialog looks fine with no missing buttons when I run Qgis itself.

Here's a test script that fires up a symbology dialog that produces the error for me, and one line after it that also produces the error when trying to load an icon. The missing icons are under the symbology layer list, with functions like new layer, lock colour, move layer up/down etc:

from qgis.core import *
from qgis.gui import *
import sys
sys.path.append("/usr/share/qgis/python/plugins")
from PyQt4.QtCore import *
from PyQt4.QtGui import *

app = QgsApplication([], True)
app.setPrefixPath("/usr", True)
app.setThemeName("default")
app.initQgis()
print QgsApplication.showSettings()

layer = QgsVectorLayer("point","point","memory")
if not layer.isValid():
    raise ValueError,"Layer not valid"

canvas=QgsMapCanvas()
dialog = QgsRendererV2PropertiesDialog(layer, QgsStyleV2.defaultStyle(), embedded=False)
dialog.exec_()

print QIcon(QgsApplication.iconPath( "symbologyAdd.svg" ) )

app.exitQgis()

Note that the settings printed out by QgsApplication.showSettings() are the same in this script and in Qgis. I'm thinking there's some app initialisation that I've not done. I can't find the resource file that has those icons in it. The icon path is the same in my script, namely: ":/images/themes/default/symbologyAdd.svg".

So, any solution?


texts:
I'm working with a GeoServer WFS in QGIS 3.4. If I open two parallel projects with the WFS layer included and successfully change some date over WFS-T in one project, I can't get the new changes in the other. Only after reopen the project or reimport the WFS layer will show the changes.

Is there an option to refresh the WFS layer after zoom and pan?

I can't find any cache options for a WFS layer.


texts: ogr2ogr conversion between from shapefile to MBtiles not working as intended

texts:
I am trying to convert a shapefile to MBTiles format using ogr2ogr and have been running this script in the OSGEO4W command line

ogr2ogr -f MBTILES -dsco MINZOOM=10 -dsco MAXZOOM target.mbtiles C:\users\myself\path_to_shp\source.shp

Yet I get this error when I run this command

ERROR 1: sqlite3_exec(PRAGMA page_size = 4096;PRAGMA synchronous = OFF;PRAGMA journal_mode = OFF;PRAGMA temp_store = MEMORY;CREATE TABLE temp(z INTEGER, x INTEGER, y INTEGER, layer TEXT, idx INTEGER, feature BLOB, geomtype INTEGER, area_or_length DOUBLE);CREATE INDEX temp_index ON temp (z, x, y, layer, idx);) failed:
ERROR 1: MVT driver failed to create target.mbtiles

This is my first attempt to use ogr2ogr and I am wondering if it is something simple I am doing wrong here.


texts: Adding custom basemap in ArcGIS.Com that can be used in Javascript API?

texts:
It's possible to add a layer as a basemap in ArcGIS.com by checking the check box use as baselayer:

But when I consume the map using:

var mapDeferred = esri.arcgis.utils.createMap("<guid>", "map", {
                mapOptions: {
                    slider: true,
                    nav: false
                }
            });

the basemap don't show up, is there any rules that define what content can be used as a basemap?


texts:
I have 2 lists (features) and I wanted to join in a graph, where the X axis would be the months and the Y axis would be the monthly data of TRMM and the other one of MODIS.
var multiYrGraph2 = ui.Chart.feature.byFeature({
  features: multyrMonthTS,
  xProperty:'month',
  yProperties: 'value'})
      .setChartType("ColumnChart").setChartType('LineChart')
      .setOptions({vAxis: {title: 'LST [deg. C]'},
                    hAxis: {title: 'Month'}});

https://code.earthengine.google.com/ea2d1f05fd2e642982057fa32d3c15bf

Actually I am trying to put the TRMM data on a Y-axis (left) and the MODIS data also on a Y-axis (right), and the X-axis would be the months.
Example:


texts: Deleting polygons that overlap specific attributes?

texts: I have a raster dataset with specific features and values. I also have a dataset that is made up of several hundred small polygons identifying specific features. How do I remove polygons that don't intersect specific values in the raster dataset?


texts:
I am working with r.stats. My analysis is about overlap Extreme Climatic Events (ECE) raster map with mammals range maps. I used to do my analysis first creating r.mask on ECE map and then applying r.stats to raster range map for single animal. I discovered I cannot use raster maps for resolution matters, so I have available only vector range maps. It means I have to invert the process and do my statistics on ECE map: so first selecting the region/mask to work on, e.g. the vector range map, and the apply r.stats on ECE raster map.
As r.mask is only usable for rasters, should I then use g.region for vector range maps? Would it work the same?


texts: Find number of points within certain radius of a set of points

texts:
I have two tables. One is a list of stores (with lat/long). The other is a list of customer addresses (with lat/long). What I want is a query that will return the number of customers within a certain radius for each store in my table. This gives me the total number of customers within 10,000 meters of ANY store, but I'm not sure how to loop it to return one row for each store with a count.

Note that I'm doing this queries using cartoDB, where the_geom is basically long/lat.

SELECT COUNT(*) as customer_count FROM customer_table
WHERE EXISTS(
    SELECT 1 FROM store_table
    WHERE ST_Distance_Sphere(store_table.the_geom, customer_table.the_geom) < 10000
)

This results in a single row :

customer_count
4009

Suggestions on how to make this work against my problem? I'm open to doing this other ways that might be more efficient (faster).

For reference, the column with store names, which would be in one column is store_identifier.store_table


texts:
I am developing an application that needs to show map tiles without internet. My initial step is to show map tiles, that part is done. Next step is to show tiles without internet. I tried to get tile images from my local directory but program stuck on the beginning and reading tiles from the given directory. What should I do? Should I serve my tiles on a local server? Or what is the best way to work at offline mode.
import QtQuick.Controls 2.4
import QtLocation 5.8
import QtPositioning 5.2
import QtGraphicalEffects 1.0

Item {
    height: 1080
    width: 1920
    Rectangle {
        anchors.fill: parent
        Plugin {
            id: plugin
            name: "osm"
            PluginParameter{
                name: "osm.mapping.offline.directory"
                value: "MAP_TILE_PATH"
            }
        }

        Map {
            id: maps
            anchors.fill: parent
            plugin: plugin
            gesture.enabled: true
            gesture.acceptedGestures: MapGestureArea.PinchGesture | MapGestureArea.PanGesture
                                      | MapGestureArea.FlickGesture
            gesture.flickDeceleration: 3000
            focus: true
            zoomLevel: 10
            center: QtPositioning.coordinate(41,30)
            minimumZoomLevel: 1
            copyrightsVisible: false
            onCenterChanged: {
                console.log(maps.center)
            }
        }
    }
}


texts: Kicking someone off a shapefile in ArcMap

texts:
I have ArcMap open, and someone else has the same shapefile open as me. How can I kick them off so that I can edit? I get the following error message:
"Unable to save edits, cannot acquire a lock"


texts: Remove empty polygons from shapefile in R

texts:
I have a shapefile witch contains several empty polygons. I´m extracting values from raster file based in some polygons of this shapefile. But my loop stops at every empty polygon.
Simple solution would be to remove NAs with na.omit(shape@data) but this deletes some xy values as well.

Here is my code. The loop stops once it finds a empty polygon.

all.mun <- NULL
  for (j in mun.codes)
  {
    mun       <- mun.br[mun.br$CODIGO_MUN == j, ] # isolate municipality
    suit.mun  <- extract(suit, mun)               # extract suits from the municipality
    suit.mun  <- do.call("c", suit.mun)           # convert to data frame no matter how many lists the function extract::raster returns
    suit.mean <- mean(suit.mun)                   # mean suitabilitie of the site
    all.mun   <- rbind(all.mun, suit.mean)        # saving result by observation, i.e municipalitie
  }

mun.condes is a vector with a list of 6 digit municipalities codes.
mun.br is the shapefile of the municipalities

> names(mun.br)
 [1] "ID"         "NOMEUF"     "REGIAO"     "CODMESO"    "NOMEMESO"   "CODMICRO"
 [7] "NOMEMICRO"  "NOME_1"     "NOME_2"     "CODIGO_UF"  "CODIGO_MUN" "AREA_1997"

Here is the error message.

Error in y[i, ] :
  cannot get a slot ("Polygons") from an object of type "NULL"


texts: Finding polygon within polygon using ArcGIS for Desktop?

texts:
How do I find two parks (shapefile) that are within a census tract ?

If I'm not mistaken you have to use selection by location.. but beyond that quite sure what to do after.


texts: How do I essentially crop an ArcGISTileMapServiceLayer as an object while maintaining the integrity of the basemap

texts:
I am using the up-to-date flex api and esri api for flex. I know how to use flex to obtain a snapshot of a map, but ultimately what I want to do is retrieve a scaled up image of a part of the basemap. Are there any tools in the flex api that I could use to modify the extent and scale of a map object and store the result maybe as a clone? I've tried changing the extent, but it ultimately changes the basemap which is not what I want to do. I simply want to select a part of the basemap using extent coordinates and then retrieve only what is within those coordinates.

I want to be able to get this image and I understand that the tiles will need to wait to be cached. So I can probably write a conditional check for that.


texts: How to split polyline using points WITHOUT splitting at intersections?

texts: I am trying to use points to split a shapefile of lines that have been dissolved.The "SplitLinesAtPoint" tool splits at intersections in addition to the points, WHICH I DO NOT WANT. I only want to split the lines at the points and remain connected at the intersections.


texts: Aggregating DateTime and GPS latitude/longitude

texts:
I have a dataframe with multiple columns including SerialNumber, DateTime, GPS latitude and GPS longitude. Background: Tractors marked with Serial numbers are working(driving) on fields, sending data about GPS latitude and longitude every few seconds. Since the data is coming in every few seconds the dataset is huge (few million rows).

EDIT:
I want to group by my data on SerialNumber and DateTime. Furthermore i want to aggregate DateTime on lets say 5min, and other features like TotalWorkingHours and AverageFuelConsumption on sum/mean/max etc. The problem comes with aggregating GPS data. From GPSLatitude and  GPS Longitude i  created geometry column, which is GeoSeries, and i want to use geometry.apply(lambda x: x.centroid) for aggregating geometry after group by.

This is my code so far but it is not working:

geometry = [Point(xy) for xy in zip(trac_df.GpsLongitude, trac_df.GpsLatitude)]
crs = {'init': 'epsg:4326'}
gd_trac_df = gpd.GeoDataFrame(trac_df, crs=crs, geometry=geometry)

trac_agg = gd_trac_df.groupby(['SerialNumber', pd.Grouper(key='DateTime', freq='5min')]).agg({
             'geometry' :  'apply(lambda x: x.centroid)', #It obviously fails here
             #'GPSLongitude' : 'median',
             #'GPSLatitude' : 'medain',
             'TotalWorkingHours' : 'max',
             'Engine_rpm' : 'mean',
             'EngineLoad' : 'mean',
             'FuelConsumption_l_h' : 'mean',
             'SpeedGearbox_km_h' : 'mean',
             'SpeedRadar_km_h' : 'mean',
             'TempCoolant_C' : 'mean',
             'PtoFront_rpm' : 'mean',
             'PtoRear_rpm' : 'mean',
             'GearShift' : lambda x: ','.join(x.astype(str)),
             'TempAmbient_C' : 'mean',
             'ParkingBreakStatus' : 'min',
             'DifferentialLockStatus' : 'max',
             'AllWheelDriveStatus' : lambda x: ','.join(x.astype(str)),
             'CreeperStatus' : lambda x: ','.join(x.astype(str))}).reindex(['TotalWorkingHours', 'geometry',
                      'Engine_rpm',  'EngineLoad', 'FuelConsumption_l_h', 'SpeedGearbox_km_h',
                        'SpeedRadar_km_h',  'TempCoolant_C', 'PtoFront_rpm',  'PtoRear_rpm',
                         'GearShift', 'TempAmbient_C', 'ParkingBreakStatus', 'DifferentialLockStatus',
                         'AllWheelDriveStatus', 'CreeperStatus', 'geometry'], axis=1)

Any suggestions how to deal with aggregating geometry GeoSeries within group by?

EDIT #2:
After using the code, i got the following message:

'SeriesGroupBy' object has no attribute 'apply(lambda x: x.centroid)'


texts:
I have produced an extensive set of annotation in ArcMap 10. I'm working with someone who basically has the need of a spatial database but only understands printed maps, so the text has grown to the point where I can't fit any more on the page. I'd like to export the annotation to a spreadsheet(no problem there) and replace the text shown on the page with an ID#. I'll then be left with a map with a bunch of small callouts which refer to a spreadsheet that my client knows how to work with. This annotation is not feature linked.

How do I change the display of my annotation feature class from the field titled 'TextString' to 'OBJECTID' or any other field?


texts:
In Arcmap, is there a plugin/addon available that will list the layers that have been edited when the "Save Edits" dialog pops up? Sometimes I get nervous that I might have deleted something from a layer that I didn't intend to edit. If so, I don't want to save my edits. Before clicking 'Save Edits', it would be helpful to review a list of layers that were edited.

I know I can toggle layer selectability, however, I have made mistakes in the past and I thought this would be an extra layer of defense against mistakes.


texts: How to secure WMS/WFS request from Geoserver?

texts:
Currently anyone can go to inspect-> network and get my rest url and access my data from geoserver, I want to secure my connection (probably by adding proxy server) so that only authenticated persons can get data when they run the url.
I have approximately 10000 users.