texts: PostGIS/GeoDjango: given a list of locations, how to find closest point of interest in a different table?

texts:
in my database on PostgreSQL/PostGIS I've two tables: a "location" table which main geo field is a MultiPolygon and a "busstop" table which main geo field is a Point.

Given a base point, I would like both to filter all the location within a certain distance from that point (I already know how to do that) and find the closest bus stop with its distance to that point.

Basically I would like to have back a recordset like this:

--------------------------------------
| location_1 | bus_stop_234 | 1.1 Km |
--------------------------------------
| location_2 | bus_stop_106 | 1.2 Km |
--------------------------------------
| location_3 | bus_stop_827 | 1.6 Km |
--------------------------------------
| location_4 | bus_stop_365 | 1.8 Km |
--------------------------------------

(Where location comes from location table and bus_stop from busstop table, while distance is calculated between multipolygon field of location and point field of bus stop.)

Even better if I can order the results by distance field, like in the example.

I need to do this using a GeoDjango query, but even a simple SQL query could be good to let me understand how to perform this kind of query.

Please note that there aren't any FK relations between the two tables.

Update: I still don't know how to do this with Django ORM and it's still not the exact query I want, but at least I was able to write a SQL query that returns the locations with the bus stops within 200 meters. I would want just the first (with closest bus stop) location/busstop for each location:

SELECT "api_location"."id", "api_location"."name", "api_location"."point", "api_location"."geom", "api_busstop"."amic_code",
ST_Distance("api_location"."geom", "api_busstop"."point") as distance
FROM "api_location", "api_busstop"
WHERE ST_DWithin("api_location"."geom", ST_GeomFromEWKB('\001\001\000\000 \346\020\000\000\210\215\325\247\202\367\001\300q\216\350\015\272\301J@'::bytea), 2000.0)
    AND ST_DWithin("api_location"."geom", "api_busstop"."point", 200)
ORDER BY id, distance


texts: R - Stack rasters with different dimensions and extent

texts:
I am trying to stack multiple rasters together in one RasterStack object but am having problems with aligning the extents and dimensions/resolution so that it can work. All rasters have the same CRS and projections.
The four raster files I am using can be downloaded from this dropbox link:
https://www.dropbox.com/sh/1fafynlmoveqr2z/AACdQlHMzwRgfdnKKxP8mcm_a?dl=0
I have tried using the following, but I am not sure if this is the best/correct method:
    library(raster)
    library(sp)
    library(rgeos)
    library(rgdal)
    library(dismo)

    # Load country administrative boundaries:
    data("wrld_simpl")

    r1 <- raster("./NDVI_metu_saudi_arabia.2021001_NDVI.MOD44CQ.250m.tif")
    r2 <- raster("./slope_1KMmd_SRTM.tif")
    r3 <- raster("./wc2.1_30s_bio_1.tif")

    # Mask raster with country border shapefile
    mskd_rstr <- mask(
      r1,
      wrld_simpl[177,]  # country border
    )

    # Trim raster to country boundary. Make all pixels outside border == NA.
    trmmd_rstr <- trim(mskd_rstr)

    # Check result.
    plot(trmmd_rstr)

    # Write to disk.
    writeRaster(
      trmmd_rstr,
      file = paste0(
        "r1",
        "_SAU_trmmd_rstr.tif"
      )
    )

# Repeat for r2, r3, r4 so that all rasters are the same and can be stacked.

This method takes extremely long (especially when using trim()). Is there a better alternative?
Ultimately, I want to be able to stack them together without getting errors about the extent or dimensions. For example by using:
predictors <- stack(r1, r2, r3, r3)


texts:
Can AIS marine traffic live vessel data be integrated into QGIS?

I understand this could be possible using an OpenCPN plugin and Grass.


texts: St_Difference between raster and vector using PostGIS

texts:
I want to calculate the difference between a raster and vector layer,

I used this querie but it's so slow and I don't know if it's right or not.

I'm not using st_difference because it only support vector layers

SELECT ST_Clip(a.raster,b.geom)
FROM my_raster a, my_shape b
WHERE not ST_Intersects(a.raster, b.geom)

I want to erase raster region wish is intersected with the vector layer

EDIT 1

Queries listed bellow work but they take a lot of time

To rasterise vector layer

create table rasterise_vector
select ST_AsRaster(st_transform(geom,4326), (select st_union(raster) from my_raster) , 2BUI') as raster
from shp

To calculate difference between rasters

create table diff_rasters as
SELECT ST_MapAlgebra(st_union(a.raster), st_union(b.raster), 'rast1 - rast2', '32BF', 'FIRST', NULL, 'rast1', NULL)
from my_raster a, rasterise_vector b

Table size :

Raster size Rows (estimation) = 1046 tiles

Vector size Rows (estimation) = 9854 items

could you please help me solve this problem?


texts:
I have a numpy array that contains floats values.
Once converted to raster, the values are all integers.
How can i keep the float format?

The procedure I use is:

geotransform = [xmin, xsize, 0, ymax, 0, -ysize]
Raster = np.array(values_array)
driver = gdal.GetDriverByName('GTiff')
outputRaster = 'output/path/to/raster.tif'
dst_ds = driver.Create(outputRaster, cols, rows,  1, gdal.GDT_Int16)
band = dst_ds.GetRasterBand(1)
band.WriteArray( Raster )
band.SetNoDataValue(-9999)
dst_ds.SetGeoTransform(geotransform)
epsg = MaskGrid.crs().postgisSrid()
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dst_ds.SetProjection( srs.ExportToWkt() )
dst_ds = None

What's wrong?


texts: My geometry appears far away from the real

texts:
I'm kind of new on PostGIS. I have the coordinates (UTM) in a table, and I need to update the geometry field based on that coordinates.

The SRID in the table is 32717. I tried this sentence
UPDATE public.uems
    SET coordenada_uem=ST_Multi(ST_SetSRID ((ST_MakePoint(172332.96000 , 68435.29000)), 32717))
    WHERE gid=78;`

And the result is this when I've tried to see the geometry

The updated point is really far away from the real ubication, and even the country (Ecuador).
How can I fix that?


texts:
I am using QGIS 2.4.0 Chugiak and am trying to create an elevation profile from a raster. When I click on the terrain profile icon I get an error message (see below).
Has these problem been encountered before and does anybody know an easy solution to this problem?
An error has occured while executing Python code:

Traceback (most recent call last):   File
"C:/Users/eesje/.qgis2/python/plugins\profiletool\profileplugin.py",
line 299, in changePlotLibrary
self.wdg.addPlotWidget(self.plotlibrary)   File "C:/Users/eesje/.qgis2/python/plugins\profiletool\ui\ptdockwidget.py",
line 145, in addPlotWidget
self.plotWdg = PlottingTool().changePlotWidget("Matplotlib", self.frame_for_plot)   File
"C:/Users/eesje/.qgis2/python/plugins\profiletool\tools\plottingtool.py",
line 82, in changePlotWidget
from matplotlib.figure import Figure   File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python\qgis\utils.py", line 454, in
_import
mod = _builtin_import(name, globals, locals, fromlist, level)   File
"C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\matplotlib\figure.py",
line 34, in
import matplotlib.colorbar as cbar   File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python\qgis\utils.py", line 454, in
_import
mod = _builtin_import(name, globals, locals, fromlist, level)   File
"C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\matplotlib\colorbar.py",
line 29, in
import matplotlib.collections as collections   File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python\qgis\utils.py", line 454, in
_import
mod = _builtin_import(name, globals, locals, fromlist, level)   File
"C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\matplotlib\collections.py",
line 23, in
import matplotlib.backend_bases as backend_bases   File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python\qgis\utils.py", line 454, in
_import
mod = _builtin_import(name, globals, locals, fromlist, level)   File
"C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\matplotlib\backend_bases.py",
line 50, in
import matplotlib.textpath as textpath   File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python\qgis\utils.py", line 454, in
_import
mod = _builtin_import(name, globals, locals, fromlist, level)   File
"C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\matplotlib\textpath.py",
line 11, in
import matplotlib.font_manager as font_manager   File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python\qgis\utils.py", line 454, in
_import
mod = _builtin_import(name, globals, locals, fromlist, level)   File
"C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\matplotlib\font_manager.py",
line 1356, in
_rebuild()   File "C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\matplotlib\font_manager.py",
line 1341, in _rebuild
fontManager = FontManager()   File "C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\matplotlib\font_manager.py",
line 1004, in init
self.ttflist = createFontList(self.ttffiles)   File "C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\matplotlib\font_manager.py",
line 584, in createFontList
prop = ttfFontProperty(font)   File "C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\matplotlib\font_manager.py",
line 398, in ttfFontProperty
sfnt4 = sfnt4.decode('ascii').lower() UnicodeDecodeError: 'ascii' codec can't decode byte 0x82 in position 0: ordinal not in range(128)
Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]


texts: Quickly browsing of features for the active vector layer?

texts:
I often need to view on the canvas the (selected or not) features from the current active vector layer.

In this situation, I only need to view something on the canvas, without doing anything on the canvas itself (so, the active window is always the Attribute Table).

This is what I usually do:

Open the Attribute Table for the current active layer;
Select the first feature;
Press Ctrl+P button (i.e. the "Pan map to the selected row" button);
Minimize the Attribute Table;
View something on the canvas for the current feature;
Go to the first step until needed.

Is there a way (a method, a shortcut) for directly browsing every feature on the canvas without the needing of going backwards and forwards with the Attribute Table? I mean, something like a "next" button?

(I may use the workaround of splitting the windows of the Attribute Table and of the canvas within the monitor, or also docking the Atrribute Table, but I would like to work on the full QGIS window).

EDIT I saw that my question was casted as a duplicate of this question. The two questions are very similiar, but I think it should be left as open for two reasons:

the related duplicated question doesn't have an accepted answer;
my question is slightly different because I'm looking for the possibility of displaying only the selected features (and the @Joseph's answer, that I marked as accepted, answers well to my request).


texts: Converting *.xlsx datasheet with point feature data into POINT feature class using ArcPy?

texts:
I am trying to create a script tool to automate this process, eventually it will be used for batch processing but at the moment I can't even get test data to work. This is pretty basic at the moment, and I know it can be done through ArcGIS, but I would like to develop a python script for processing. My code is as below: It completes without error, but I am expecting the creation of a feature layer (even better would be a feature class with POINT geometry). I can't find or access the output of the SaveToLayerFile step.

#Import libraries
import arcpy
import os
import sys
arcpy.env.overwriteOutput = True

#Set workspace
home = arcpy.env.workspace
arcpy.env.workspace = arcpy.GetParameterAsText(0)
mxd = arcpy.mapping.MapDocument("Current")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]

#Create spatial reference object
prjFile = arcpy.GetParameterAsText(2)
spatialRef = arcpy.SpatialReference(prjFile)
arcpy.AddMessage('Setting Spatial Reference')

#Set Excel Datasheet for Migration
TEST = arcpy.GetParameterAsText(1)

try:
# Process: Table to Table Conversion with destination variable
    FarmerTable = 'Farmers1'
    arcpy.AddMessage('Migrating Data...')

    arcpy.TableToTable_conversion (TEST, home, FarmerTable)
except:
    arcpy.GetMessages()

# Set local variables
in_layer = "FarmerLyr"
out_layer = "FarmerLyr.lyr"
out_layer0 = "FarmerLyr"

try:
#Make XY Events of each Farmer
    arcpy.AddMessage('Creating Features...')

    result = arcpy.MakeXYEventLayer_management(FarmerTable, "X", "Y", out_layer0, spatialRef, "")
    result.getOutput(0) = out_layer0
    arcpy.AddMessage(arcpy.GetCount_management(out_layer0))
    arcpy.AddMessage('Successfully Created Features')

    # Execute SaveToLayerFile
    arcpy.SaveToLayerFile_management(out_layer0, out_layer, "ABSOLUTE")
except:
    print arcpy.GetMessages()


texts: Adding multiple feature classes (gdb) to map file (aprx) by using addDatafromPath in ArcPy?

texts:
I have some problems adding multiple feature classes (not only a single one) from a geodatabase to a map file in an ArcGIS Pro project and I'm wondering if someone knows of a solution to this.
To be able to add more than just one layer from the gdb I've tried to use ListFeatureClasses in order to select all feature classes before using AddDatafromPath. But this doesn't seem to do the trick either.

Here's the script. I've tried several other solutions as well, and while some of them did not return any errors they still did not add the data to the map file.

arcpy.env.workspace = "path to gdb"
datasets = arcpy.ListDatasets(feature_type='feature')
datasets = [''] + datasets if datasets is not None else []
rootPath_aprx = "path to aprx"
aprx = arcpy.mp.ArcGISProject(rootPath_aprx+ 'Test.aprx')
aprxMap = aprx.listMaps("Map")[0]
for ds in datasets:
    for fc in arcpy.ListFeatureClasses(feature_dataset=ds):
        path = os.path.join(arcpy.env.workspace, ds, fc)
        print(path)
        try:
            aprxMap.addDataFromPath(test)
            print("yeah.")
        except:
            print("gah.")


texts: Calculating head drop along river automatically in ArcGIS for Desktop?

texts:
I am calculating the hydropower potential of a large river basin. I have already had the hydrological modeling of the large (53000 km2) river basin and also I have DEM of the basin  and river network. I am trying to calculate the cascade of hydropower potential along the river network. I want to start at any downstream point of a river (for hydropower station) and move upstream along the river in such a way that the head drop between the initially start point and final point would be greater than or equal to 10 m. When the head drop is more than or equal 10 m, I would like to move 1000 m upstream for next hydropower station assuming this as the second point for hydropower station. From this second start point I want to move upstream such that the head drop would be greater than or equal to 10 m. Likewise I would like to repeat along the whole river network till the end of river.

I could not work out it. Please help me how can I do it automatically and manually as well. As my river basin is very large (53000 square kilometer) , I want to do it automatically in ArcGIS.


texts: Creating bounding box in ArcGIS API for JavaScript?

texts:
I am trying to create a bounding box for my map in ArcGIS 3.28. I want the bounding box to stay in a given area based on the coordinates given. I found a class called Extent in the API Reference that can set a bounding box, but it does not look like there's a way to implement the bonding box. This is what I have at the moment:`

var map;

  require([
    "esri/map",
    "esri/dijit/InfoWindowLite",
    "esri/InfoTemplate",
    "esri/config",
    "esri/layers/FeatureLayer",
    "dojo/dom-construct",
    "esri/geometry/Extent",
    "esri/dijit/Search",
    "esri/layers/ArcGISTiledMapServiceLayer",
    "esri/SpatialReference",
    "dojo/domReady!"
  ],

  function(
    Map,
    InfoWindowLite,
    InfoTemplate,
    esriConfig,
    FeatureLayer,
    domConstruct,
    Extent,
    Search,
    ArcGISTiledMapServiceLayer,
    SpatialReference
  ){

     var boundingBox = new Extent(-122.3737, 37.454, -121.467, 37.9067, new SpatialReference({ wkid:4326 }));

       //creating the map
     map = new Map("map", {
       extent: boundingBox,
       basemap: "streets",  //For full list of pre-defined basemaps, navigate to http://arcg.is/1JVo6Wd
       center: [-121.91, 37.65], // longitude, latitude
       zoom: 11
     });

Using this, is there a way to create a bounding box or is something I can add to it to make it work?


texts: Is there a free tool for configuring WMS configuration xml?

texts: Instead of writing WMS xml file from scratch with notepad , is there a tool that is easy to create/configure such a file (or convert a GIS workspace - MapInfo / QGIS into WMS file) ?


texts: Why are Moran's Index calculated in GeoDa and ArcGIS different?

texts:
I recently have found that the Moran's Index calculated in GeoDa and ArcGIS are different based on the same data and the same spatial matrix.
For example, this is my testing shapefile, a point shapefile with 547 records and UTM-Zone-50N projection. Here's display of the data and records.

I created 4-nearest neighbors spatial matrix and calculated the Moran's Index of variable 'lst', using GeoDa 1.6.7, 1.8.12 and ArcGIS 10.2. Results should be absolutely the same, but I got 0.672521, 0.674407, and 0.67824, respectively. Although the difference is very small, we cannot ignore it.  And And I am sure the data is fine, no missing/undefined value.

Does anybody find similar problem? Or is there anyone who knows about this problem?

I shared the test data in my google drive (https://drive.google.com/open?id=0BziPQVV2jY0PU2V6Ulo2WE1zeDQ).

This is my python script to calculate Moran's Index.

import arcpy
from arcpy import env
import os
env.overwriteOutput=True
pointfile="e:/temp.shp"
pointswm ="e:/temp.swm"
try:
    #4-nearest neighbors spatial matrix
    swm = arcpy.GenerateSpatialWeightsMatrix_stats(pointfile,"num",pointswm,
                                                "K_NEAREST_NEIGHBORS",
                                                "EUCLIDEAN", "1", "", "4",
                                                "ROW_STANDARDIZATION",
                                                "", "", "", "")
    #calculate Moran Index
    mo = arcpy.SpatialAutocorrelation_stats(pointfile,"lst","NO_REPORT",
                                        "GET_SPATIAL_WEIGHTS_FROM_FILE",
                                        "EUCLIDEAN_DISTANCE", "NONE","",
                                        pointswm)
    print mo.getOutput(0)
except arcpy.ExecuteError:
    print arcpy.GetMessages()


texts:
I have been trying to create a WKT entry to use in gdalwarp etc. I need an affine transform from coordinates in MGA94 Zone 54 to our local mine grid.  All efforts have so far failed, including using the bits and pieces I have found here.

I have a MapInfow.prj file entry that works beautifully but I need to convert our imagery from MGA to mine grid to supply to mining consultants. This entry is below with the affine parameters highlighted:

"SRG/AISTR", 3008, 116,7, 141, 0, 0.9996, 500000, 10000000, 7,0.999770876952, -0.031210618094, -193025.578125, 0.031224025413, 0.999768821523, -7589075.9140625, -10000, -10000, 500000, 500000

The standard scale, rotate parameters as supplied by our surveyors are here:

Rotation: 358.2111111 degrees
Scale: 0.999741893
Local East: 16000
Local North: 20000
MGAz54 E: 446231.629
MGAz54 N: 7596899.020

My latest attempt using the affine parameters from MapInfo are here:

FITTEDCS["SRG",
   PARAM_MT["Affine",
       PARAMETER["num_row, 3],
       PARAMETER["num_row, 3],
       PARAMETER["elt_0_0", -193025.578125],
       PARAMETER["elt_0_1", -7589075.9140625],
       PARAMETER["elt_1_0", 0.999770876952],
       PARAMETER["elt_1_1", 0.031224025413],
       PARAMETER["elt_2_0", -0.031210618094],
       PARAMETER["elt_2_1", 0.999768821523],
       PARAMETER["Central Meridian", 140.4800753],
       PARAMETER["Latitude_Of_Origin", -21.73061814]],
   PROJCS["GDA_1994_MGA_Zone_54",
       GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
       PRIMEM["Greenwich",0.0],
       UNIT["Degree",0.0174532925199433]],
   PROJECTION["Transverse_Mercator"],
   PARAMETER["False_Easting",500000.0],
   PARAMETER["False_Northing",10000000.0],
   PARAMETER["Central_Meridian",141.0],
   PARAMETER["Scale_Factor",0.9996],
   PARAMETER["Latitude_Of_Origin",0.0],
   UNIT["Meter",1.0],AUTHORITY["EPSG",28354]]

As requested by Andre here are some control points:

Local East    Local North MGAz54 East MGAz54 North
15896 27965   446376.2852 7604861.329
12458 28452   442956.0464 7605455.262
35698 27412   466146.2673 7603690.743
15698 29654   446231.1445 7606555.249
16547 23456   446886.0792 7600335.373


texts: I am trying to get data from table A to table B. But it is important that domains in table B would remain. So I tried to use Append, but it doesn't work with tables. Any ideas?


texts:
I install QGIS via OSGeo4W 32 bit on windows 7 32 bit and i want to take PYQGIS API in standalone python script.

First I have python 3.5.2 and I have before install  sip via pip and pyqt4 via gohlke booth successful.]]

After install QGIS via OSGeo4W and I create a qgis.pth file in python site-packages with the paths.

sys.path
['', 'C:\\Users\\usename\\AppData\\Local\\Programs\\Python\\Python35-32\\Lib\\idlelib', 'C:\\OSGeo4W\\apps\\qgis\\python', 'C:\\Users\\usename\\AppData\\Local\\Programs\\Python\\Python35-32\\python35.zip', 'C:\\Users\\usename\\AppData\\Local\\Programs\\Python\\Python35-32\\DLLs', 'C:\\Users\\usename\\AppData\\Local\\Programs\\Python\\Python35-32\\lib', 'C:\\Users\\usename\\AppData\\Local\\Programs\\Python\\Python35-32', 'C:\\Users\\usename\\AppData\\Local\\Programs\\Python\\Python35-32\\lib\\site-packages', 'C:\\OSGeo4W\\bin', 'C:\\OSGeo4W\\lib', 'C:\\OSGeo4W\\apps', 'C:\\OSGeo4W\\apps\\qgis\\bin', 'C:\\OSGeo4W\\apps\\qgis\\lib', 'C:\\OSGeo4W\\apps\\qgis\\qtplugins', 'C:\\OSGeo4W\\apps\\qgis\\python\\qgis', 'C:\\OSGeo4W\\apps\\qgis\\resources', 'C:\\OSGeo4W\\apps\\qgis\\plugins', 'C:\\OSGeo4W\\apps\\qgis\\grass\\bin', 'C:\\OSGeo4W\\apps\\qgis\\grass\\modules', 'C:\\OSGeo4W\\apps\\qgis\\grass\\scripts', 'C:\\OSGeo4W\\apps\\saga', 'C:\\OSGeo4W\\apps\\saga\\dll', 'C:\\OSGeo4W\\apps\\saga\\modules', 'C:\\OSGeo4W\\apps\\gdal', 'C:\\OSGeo4W\\apps\\gdal-dev\\bin', 'C:\\OSGeo4W\\apps\\grass\\grass-7.3.svn\\bin', 'C:\\OSGeo4W\\apps\\grass\\grass-7.3.svn\\etc\\python', 'C:\\OSGeo4W\\apps\\grass\\grass-7.3.svn\\etc\\python\\grass', 'C:\\OSGeo4W\\apps\\grass\\grass-7.3.svn\\lib', 'C:\\OSGeo4W\\apps\\grass\\grass-7.3.svn\\scripts', 'C:\\OSGeo4W\\apps\\Python25', 'C:\\OSGeo4W\\apps\\Python27']

but if i try to import qgis.core show me that error message :

Traceback (most recent call last):
  File "<pyshell#3>", line 1, in <module>
    import qgis.core
  File "C:\OSGeo4W\apps\qgis\python\qgis\__init__.py", line 26, in <module>
    from qgis.PyQt import QtCore
  File "C:\OSGeo4W\apps\qgis\python\qgis\PyQt\QtCore.py", line 60, in <module>
    QPyNullVariant.__nonzero__ = MethodType(__nonzero__, None, QPyNullVariant)
TypeError: method expected 2 arguments, got 3

Any idea why ?


texts: Why can't I transpose a raster across the pole with rasterio or gdal, and how could it be done?

texts:
I have a raster in e.g epsg:4326, it covers all world, from 180E to 180W, from 90S to 90N. I want to re-project a part of it to epsg:3031 and also change the grid size and resolution.

Here is a minimal, and working example:

from rasterio.transform import from_bounds, from_origin
from rasterio.warp import reproject, Resampling
import rasterio
import matplotlib.pyplot as plt
import numpy as np

km = 1000

my_raster = rasterio.open("my_raster.tif") #Global raster
west, south, east, north = -180, -90, 180, 90
src_crs = my_raster.crs
src_shape = src_height, src_width = my_raster.shape
src_transform = from_bounds(west, south, east, north, src_width, src_height)
source = my_raster.read(1)

dst_crs = {'init': 'EPSG:3031'} # WGS 84 / Antarctic Polar Stereographic
print(src_crs['init'])

dst_transform = from_origin(0*km, 3000*km, 50*km, 50*km)
print(dst_transform)

my_array = np.zeros((120, 120))
my_array[:] = np.nan

reproject(source,
          my_array,
          src_transform=src_transform,
          src_crs=src_crs,
          dst_transform=dst_transform,
          dst_crs=dst_crs,
          resampling=Resampling.bilinear)

plt.imshow(my_array)

plt.show()

This raster only contains the eastern hemisphere, but I want to include all Antarctica.

When I change the transform, to e.g. dst_transform = from_origin(-1000*km, 3000*km, 50*km, 50*km) the code freezes. I've tried a number of other transformations, this is the only case that doesn't work, but I need to have my raster in Antarctic stereographic.

How can I redefine my transformation matrix so that my reprojected raster contains also Western Antarctica?

I prefer to use rasterio, as it works good for the rest of my project, but any suggestion is welcome.

UPDATE:
I tried gdal warp with Python wrapper, but also no success. Now, I can define the right bounds, but I got a hole around the South Pole instead.

A minimal (not) working example:

from osgeo import gdal
import matplotlib.pyplot as plt
import numpy as np

km = 1000

my_array = gdal.Warp('', 'my_raster.tif',
            format = 'MEM',
            width = 120, height = 120,
            dstSRS='EPSG:3031',
            outputBounds= (-2000*km, -2000*km, 2000*km, 2000*km),
            resampleAlg = gdal.GRIORA_NearestNeighbour).ReadAsArray()

ds = None

plt.imshow(my_array[1,:,:])
plt.show()


texts: Finding the average of a raster across a given region, in R, for many points in time

texts:
I have a time-series of rasters (daily maps of global temperature for the past year), and a shape-file outlining state boundaries.
I'd like to find the average temperature in each state on each day of the previous year, for a result like:
Alabama 01-01-2019 : 5.2 degrees
Alabama 01-02-2019 : 3.6 degrees
...
Wyoming 12-30-2019 : -3.1 degrees
Wyoming 12-31-2019 : -2.0 degrees

How would I do this in R? (Other project requirements prevent me from using Python or QGIS.)
The format I've tried is:
shape = st_read("./states_shp")
temp_raster = brick("temperature_data.nc")
mean_temps = exact_extract(temp_raster, shape, 'mean')

This successfully masks my raster data using the shape file, but the 'mean' operation performs a temporal mean. Is there a way I can modify this to get the spatial mean across each state while keeping the time dimension?
Edit: following a suggestion below, here's the data. The shapefile is from the US Census (https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html):
> library(sf)
> shape = st_read("cb_2018_us_state_500k/cb_2018_us_state_500k.shp")
> head(shape)
Simple feature collection with 6 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -103.0026 ymin: 28.92861 xmax: -75.24227 ymax: 40.6388
geographic CRS: NAD83
  STATEFP  STATENS    AFFGEOID GEOID STUSPS           NAME LSAD        ALAND
1      28 01779790 0400000US28    28     MS    Mississippi   00 121533519481
2      37 01027616 0400000US37    37     NC North Carolina   00 125923656064
3      40 01102857 0400000US40    40     OK       Oklahoma   00 177662925723
   AWATER                       geometry
1  3926919758 MULTIPOLYGON (((-88.50297 3...
2 13466071395 MULTIPOLYGON (((-75.72681 3...
3  3374587997 MULTIPOLYGON (((-103.0026 3...


texts:
Is it possible to express a 4D range query in Postgis? I would like to know what 4d points lie within a hypercube.
This is my table:

    CREATE TABLE myTable
    (
      point geometry,
      attrib double precision,
      CONSTRAINT ge_idx PRIMARY KEY (point)
    );

This is how I inserted the data points:

    insert into myTable (point) values (ST_MakePoint(-71.10, 42.31, 12.43, 54.3));


texts:
I wonder if someone can clarify as I'm confused.

I need to load planet OSM data into PostGIS so I downloaded a small data set of a UK County and followed the instructions here  loading the Google Mercator projection 900931.sql first.

looking at the data in GeoServer the co-ordinate readout on the location bar appears to be in metric East / West and not Longitude/Latitude yet the osm wiki states World Geodetic System 1984 is the coordinate system of OSM data.

What am I doing wrong?

Also is there a best practice for displaying OSM as a base layer?

Where can I get a Point.sld from?


texts:
I'm trying to use gdal_proximity to create a distance raster. When I try to run it independently, it works fine. However, when I try to use it within modeler, it encounters an error. I've attached two images of the dialogue boxes for the Gdal_Proximity, the first one is the dialogue when using it inside the modeler. The only difference I see is the Band number part.

This is the error message I get when running in the modeler:

  Prepare algorithm: gdal:proximity_1
  Running Rast_proximity [6/6]
  Input Parameters:
  { BAND: 1, DATA_TYPE: 5, INPUT: '', MAX_DISTANCE: 0, NODATA: 0, OPTIONS: '', OUTPUT: >'C:/Users/Reymar.Victoria/AppData/Local/Temp/processing_7b1696077f6c49b0aa61746e14dc2547/6e2edaf05ea14791966e7f34ab933e2b/gdal_proximity_1_out_rast3_proximity.tif', REPLACE: 0, UNITS: 0, VALUES: '100' }
  Traceback (most recent call last):
  File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\gdal\GdalAlgorithm.py", line 126, in processAlgorithm
  commands = self.getConsoleCommands(parameters, context, feedback, executing=True)
  File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\gdal\proximity.py", line 141, in getConsoleCommands
  raise QgsProcessingException(self.invalidRasterError(parameters, self.INPUT))
  _core.QgsProcessingException: Could not load source layer for INPUT: invalid value


texts: Programmatically grouping multiple layers into separate groups using QGIS?

texts:
I'm trying to group a large number of layers together based on a common portion of name string. They are a mix of polygons and points. I would like to create individual groups by using part of the layer's name. Anyone have any ideas? I was looking to build a plugin but I am unable to (due to my coding abilities).

I attached a picture to show what I'm looking for as the end result. I have projects with a large amount of information that needs to be grouped.


texts: ARCGIS Python script (needs a raw input line to let the user choose the output file name)

texts:
I just need a quick fix here.
Having trouble finding a way to allow the user to input the new output file name (the one that is created after the buffer function completes). Any ideas?
I am also not sure if I created the variables correctly.

Using the latest (ArcGIS 10.1 SP1)

python script to change point layer to a 50 m buffer layer
change the output file name to whatever so it can make a new one

import arcpy
arcpy.env.workspace="C:/Users/oxana/Desktop/David"
arcpy.Buffer_analysis("builteditptFT", "c:/users/oxana/desktop/david/PDpyth1","50   meters","FULL","ROUND","LIST")

note that this only buffers from the point layer that was converted from
the centroids of the newly created buildings
it is prone to generalizing around each point because the buffer radiates
50 meteres from the center of each building (not the outsides)


texts:
I have been working with the algorithm gdalogr:warpreprojectand have failed several times in the parameter in the title. Now, the problem is solved. It was about the EXTRA parameter. I read 3 books, and I know the explanation, but what I do not understand is the fact that it accepts an empty string. Furthermore, it ONLY accepts that empty string. My code is:

reproj = processing.runalg("gdalogr:warpreproject", raster, "EPSG: 25830", "EPSG:4326", '0', 0, 0, 0, 0, 2, 1, 1, False, 1, False, "", "G:\\pyQGIS\\MIRAR_COMP\\nuevo.tif")

Can anybody explain how does this parameter work?


texts:
I have a personal GDB table that has a non-unique SLID field and a date field.  Basically there are multiple records with the same SLID numbers but different dates.  I want to select only the most recent records for each SLID number.  I have tried this query statement in the Select by Attributes window:

    [Date_] IN (SELECT MAX([Date_]) FROM table_name GROUP BY [SLID])

This hasn't worked, because for some of the SLIDs, it selects multiple instances of the same SLID value, even if the dates are different.  For instance, for SLID 27, the query selects both the instances from 5/28/2013 and from 12/12/2008.  What am I doing wrong?

The SLID field is short integer format, and the Date_ field is date format.


texts: Using all dimensions with postgis functions

texts:
I have points as geometry format in table that I use for some analysis. All points are saved in four dimensions eg. x, y, z and m - lon, lat, alt and time. The issue is that for the analysis are used functions like ST_GeomFromEWKT and ST_MakeLine, which included only two dimensions form the four present.

My specific questions is can I have the four dimensions using these functions?

I tried simple example

SELECT ST_AsText(ST_MakeLine(ST_MakePoint(1, 1, 0, 1339590840), ST_MakePoint(2, 2, 0, 1339590840)));

Which returns "LINESTRING(1 1,2 2)", without z and time dimensions.

Thanks in advance!


texts:
I have DEM hillshade layers which I have tried to change the base heights to floating surfaces.

I have done this with these layers before but had to change them back and forth to do some analysis.

They keep disappearing and not coming back.

Can anyone tell me how to stop this from happening?


texts: QGIS difference processing returns unknown layer

texts:
I am trying to do a difference between two vector layers in QGIS.  I have had success in one case but I get a result I don't understand with some different data. A new item shows up in the 'Layers' field but there is no polygon in the map display. I think the layer is an 'Attribute table' but I don't recognize the icon.  Below is a screen shot of the QGIS 'Layers' area - diff is the unknown output.

What is this icon?  Any ideas why the difference operation is returning this?
Thank you.


texts:
I went to refresh the share properties on a feature service this morning and when I got to our Portal site all items were missing from the contents tab as well as in all groups. We now seem to have no content at all on the site. On Portal admin page I can see that the items still exist and are showing up in the applications pointing to their rest services. I’ve looked online for a solution which entailed reindexing and restarting the Portal for ArcGIS service as well as restarting the entire server. This did not resolve the issue. To note, the count of data items in datastore and count of index do not match even after reindexing (index counts are all at 0).
Any other ideas of how I can get the items to show back up on our Portal site?


texts: Selecting a line segment with an azimuth value closest to the azimuth value of another line segment in QGIS

texts:
I have a road network dataset and a dataset that consists of multiple disconnected lines of varying lengths in QGIS.

Most of these lines are (largely) parallel to the road network, but their shape and distance between these lines and the road network can differ. Since these lines tell something about the condition of the road segments, I would like to select/extract the road segments that are parallel to these lines. So far, I have calculated the line bearing/azimuth of each road- and line segment and added this to a new column to both attributes tables. I now need to select the road with azimuth closest to the azimuth of the disconnected line (let’s say within a percentage difference of 5 or 10%).

Any thoughts on how to do this in QGIS?


texts:
I wanting to reproject the vertices of my polygon. I have been able to achieve this for the bounds but having a bit of trouble with the vertices.

           function onFeatureModified(event) {
           var bounds = event.feature.geometry.getBounds();
           var answer = "bottom: " + bounds.bottom + "\n";
           answer += "left: " + bounds.left + "\n";
           answer += "right: " + bounds.right + "\n";
           answer += "top: " + bounds.top + "\n";
           alert(answer);

           var mapProjv =  map.getProjectionObject();
           var epsg4326v = new OpenLayers.Projection("EPSG:4326");
       var verticesNative = event.features[0].geometry.getVertices();
       var verticesLatLon = verticesNative.transform( mapProjv , epsg4326v  );
       var vertices = verticesNative;
           alert(verticesNative);
           }

Below code that works:

           var vertices = event.features[0].geometry.getVertices();
           alert(vertices);


texts: Get Numpy arrays from Google Earth Engine

texts:
I would like to get the numpy arrays of some Google Earth Engine collection.

For example here my intent would be to obtain a 4d numpy array with dimensions [timestamp, height, width, bands] (or whatever order) with the Sigma0 data of a given region, defined by its bounds (here in UTM projection).

e.g.

import ee
import geopandas as gpd
from shapely.geometry import box

ee.Initialize()

crs = {'init': 'epsg:32628'}
poly = gpd.GeoSeries(box(578150.0, 1572350.0, 578350.0, 1572550.0), crs=crs)

poly_json = eval(poly.to_json())
gj_geom = poly_json['features'][0]['geometry']
geom = ee.Geometry(gj_geom, poly.crs['init'][5:])

collection = (ee.ImageCollection('COPERNICUS/S1_GRD')
              .filter(ee.Filter.date('2019-01-01', '2019-02-1'))
              .filter(ee.Filter.geometry(geom))
             )

How do I obtain the raw data arrays from the filtered collection?


texts: Publish cached map service using pre-existing cache folder

texts:
We use:
ArcGIS Server Advanced Enterprise 10.3                                                                          /
                                                                                ArcGIS Desktop Basic 10.3
both located on a Windows Server 2012 machine

We've received a directory containing a cached map (bundle files). We don't have a mxd file, we don't have the data, just the cache.
This directory is not located on the deafault C:\arcgisserver\directories\arcgiscache, and we don't intend to move it there, as there woudn't be enough space on disk C:\ . The cache is quite large 280 GB. So I've created a new cache directory in ArcGIS Server Manager that leads to the exact path of our cache.

We would like to publish this cache as a cached map service/image service on the server. I've added it to ArcMap, but when publishing it as a map service, ArcMap doesn't automatically recognize that it's a cached service. As it is shown in this video: http://video.arcgis.com/watch/964/publishing-cached-map-services-using-a-pre_dash_existing-cache-folder

It works as a dynamic service

What is this correct procedure to make this existing cache available on ArcGIS Server as a service?


texts: Create a raster map from more than one shape

texts:
I've some shapefile loaded in Arcmap (steets, buildings...) with an appropriate symbology ad I need to convert them to unique raster, to obtain a map (like topographic map) with a specified resolution (cellsize) ad scale.

There's a method, or plugin or program to do that?

Thanks!


texts:
I want show different imagery backgrounds in a QGIS atlas based on which one is better.

I have the following highlighted imagery layers

In the dataset used for the atlas I have

How do I set up the lock layers expression in the Map Properties to just turn on the background imagery that's listed in the background column? The expected input is list of map layer names separated by | characters but I would prefer to define just the layers to add (so that all the vector data that's in the map view are also retained and just the listed background is added).


texts: Using embedded SVG symbol in multiple layers

texts:
When I embed an SVG symbol in QGIS 3.x, the symbol is stored in the layer definition. If I have another layer and like to use the symbol again, I have to embed it a second time.
How can I embed an SVG symbol globally and use it in multiple layers?


texts: Shapely : cascaded_union seems to return too many polygons

texts:
I'm trying to get the union of around 4,000 polygons using Shapely's cascaded_union function. The output of this function is a multipolygon in my case, which is normal since many input polygons do not overlap.

Yet, when I draw the output, I can clearly see that some small polygons are included whereas these polygons belong to a larger polygon. To my mind, these small polygons should be excluded from the union.

Here is a sample screenshot with 3 polygons of the union returned by shapely. The blue one is the bigger one with 43 185 exterior points and the green and orange ones are the smaller ones.

I also uploaded a text file containing the coordinates of the 3 Polygons using the built-in function mapping() from the shapely.geometry module.

It can be accessed there.

The file contains a dictionary {index:{...}} with 3 indexes, each of them containing a dictionnary representing a Polygon. If you need to get back the coordinates from the file, you can use the following Python code:

with open("file_polygon.txt","r") as f:
    testl = f.readlines()
teststr = "".join(testl)
testd = json.loads(teststr)
for index in testd:
    # List of exterior points of the polygon of index "index"
    ps = testd[index]["coordinates"][0]

If you have any clue of why these polygons are returned by the cascaded_union function, I would be happy to get them :).


texts: Creating Elevation Profile from Line Shapefile of Contours using QGIS?

texts:
I have a line shapefile that contains contour lines. The attribute table of the shapefile has a field that gives the elevation of each line. I would like to draw a line across the contours and get an elevation profile. I found the Profile tool but it seems to require some sort of raster layer, and not my line shapefile layer (gives message "active layer is not a profilable layer"). It may be that I need to use some other plugin to convert my shapefile layer to the proper layer to use in the Profile tool.

Is there some sort of step-by-step tutorial on what plugins I need and how to convert a shapefile layer to the correct type of layer so that I can obtain an elevation profile by drawing a line across the contours?

Using QGIS 2.18.18


texts: how to find optimum number of college bus routes for connecting 120 bus stops, distributed over large road network

texts: I have 120 bus stop locations with road network  and number of students. I prepared  the network dataset, but i don't know what to do next to find out optimum number bus routes which required to connect the bus stops based on the number of students using ArcGIS 10.1.


texts:
I used st_join() to combine t1 and t2.
library(sf)
j <- st_join(
  x = t1,
  y = t2
)

t1 represents voting precincts boundaries and t2, school zone boundaries. In some cases, the boundaries of t1 fit perfectly in one school zone but because the borders touch two other zones, st_join lines up a precinct with three zones.
How can I make st_join stricter so that if a precinct in t1 is mostly or seemingly entirely within a zone boundary of t2, the precinct lines up with only that zone, no others?


texts:
I'm trying to figure out how to best get started on a web-mapping project, and am wondering about where I should host my maps from.

In terms of my skills, I'm quite good at desktop mapping and cartography, and I've got mediocre PHP, mySQL, HTML and CSS skills - but I want to create this really fantastic interactive web map site that will (in my dreams only, probably!) have some kind of base layer, then overlay my own township and section grid, and have possibly hundreds of thousands of sample points. Each of these points will have data attached, and when clicked, will display it's info in a pane, and then provide the ability for registered users to comment on and discuss it. I'd like to also get some kind of feedback (thumbs up/down) ability for each user comment. Maybe I'm too optimistic, but this is what I'd like to have in the end.

My reading tells me I need a stack containing something like postGIS, MapServer, and OpenLayers but I use inmotion for my current web development and learning, which I really like, so rather than waste time figuring out how to install all of these packages there (if I even can), is it better to go with an outfit like mapserverpro.com, which has this stuff set up already and specializes in GIS hosting?

The price of their hosting packages is not a problem, and I just want to get pointed in the right direction to start with.

I guess that my short question is: does anyone here have an opinion on what mapserverpro.com is like?

Even mapserving.com?

Or would a generic hosting service usually be able to handle the stack that I mentioned, and efectively be able to handle what I picture as being a large amount of GIS data?

Also, if anyone has any comments on how to most efficiently tackle a project like this, then please feel free to set me straight.


texts:
In ARCGis its possible to Clip your Data Frame so you just see the Data within the defined area. Like shown here:
https://www.cdc.gov/dhdsp/maps/pdfs/GIS-Tips-Data-Frame-Clip-Function.pdf

With this functions the loading of large WMS Basemaps Speeds up.

Is this possible in QGIS Map View?


texts: How to print out a road atlas with QGIS?

texts:
This is what I managed to do so far (first picture):

Through the web tab and quickmapservice or OpenLayers plugin I loaded maps where on I managed to add a vector grid under research tools to then go to a composer manager where I implemented the Atlas function.

I was successful in dividing the city map into 1km square parts, but when I print I end up with a map that has the correct zoom level and quality, but the roads and street/object names have shrunk and are too small to read (second picture).


texts: I've just saved a GeoTiff image via GDAL (C++ version). I've called BuildOverviews on it. How can I check if the pyramids are actually created in the image?


texts: How to make Maximum values of a stacked raster to form a new raster in R using raster package


texts: numbering of features when dissolving non-contiguous polygons

texts:
I have dissolved polygons together to make feature A. Some of the polygons are contiguous and some are not. Therefore, there is a unique row for each feature (singlepart to multipart).

In ArcGIS or SAS how are the features assigned a unique ID?

Is the largest feature given the first ID?


texts: GeoExplorer not display the associated table to a layer

texts:
GeoExplorer shows the table information associated with a layer of two ways:

Information button: clicking on the map is displayed a tab corresponding to the polygon attributes where you click.
Query Button: The table shows the data for the entire map.

In the project that we are developing, we start with a shapefile,  manipulate it in QGIS and layers are exported in shapefiles, then uploaded onto GeoServer.

A few months ago we uploaded to GeoServer a set of maps and now we are creating new maps. The old maps being displayed on the GeoExplorer show well the tab and the table, but in the new ones neither the tab nor the table are working, only the table header displays.

Apparently the generation process in both steps the layers was the same. What can be causing this problem?


texts: How can I create a 10 day composite image of the snow cover extent using QGIS?

texts:
I have 10 singleband rasters of the snow cover extent (1 raster/day from 1st of December to 10th of December). I have to create a composite image that merges all of them. I tried the r.composite tool from GRASS and the Build Virtual Raster tool, but I don't really understand the results. The raster values are from 20 to 200:
20 - sea,
21 - lake,
22 - river,
30 - cloud,
100-200 - % of the SCE (example: 154 = 54% snow cover).
One of the input rasters:


texts: Is it possible to print WMTS tile layers with desired resolution / zoom level in QGIS?

texts:
I have a project which uses WMTS and some vector layers in QGIS. I can freely zoom in and out, and tiles are refreshed with appropriate zoom level, just as expected.

Now I would like to export part of this map (a given bounding box) with the tiles at a fixed zoom level, so that I can better match the printer resolution when map is printed in the desired physical size.

How should I tell QGIS to use that given zoom level, and which file format is best suited to do what I want?


texts: Black theme ("Night Mapping") doesn't load

texts:
At work we are using QGIS 3. When we were using QGIS 2 I enjoyed working with the black theme, called Night Mapping, for healthy eyes, etc.
With QGIS 3, when I click on the menu of the themes and pick Night Mapping, nothing happens. Am I the only one ? I want it back.


texts: Forming a single TIFF image of a Sentinel 2 tile

texts:
I downloaded the data using the following command from Amazon S3 bucket :
aws s3 cp s3://sentinel-s2-l2a/tiles/51/H/VC/2020/10/19/0/ aus1 --recursive --request-payer.

Then translated the 13 bands using :
gdal_translate ...R10m/B02.jp2 ...B02.tif
gdal_translate ...R10m/B03.jp2 ...B03.tif
gdal_translate ...R10m/B04.jp2 ...B04.tif
gdal_translate ...R10m/B08.jp2 ...B08.tif
gdal_translate ...R20m/B05.jp2 ...B05.tif
gdal_translate ...R20m/B06.jp2 ...B06.tif
gdal_translate ...R20m/B07.jp2 ...B07.tif
gdal_translate ...R20m/B8A.jp2 ...B8A.tif
gdal_translate ...R20m/B11.jp2 ...B11.tif
gdal_translate ...R20m/B12.jp2 ...B12.tif
gdal_translate ...R60m/B01.jp2 ...B01.tif
gdal_translate ...R60m/B09.jp2 ...B09.tif

gdal_merge.py -init 255 -o output.tif B01.tif B02.tif B03.tif B04.tif B05.tif B06.tif B07.tif B08.tif B8A.tif B09.tif B11.tif B12.tif

What I want to do after this is supervised classification of these images, classifying them as either land or sea or both.
My question is the TIFF files that I am not sure that the TIFF files I have produced combining all bands is a correct representation for the tile or not.


texts:
I am building an application primarily with the ESRI API, but would like to load a mapbox tilelayer to use as the basemap as they are just so much more pleasing cartographically. I dont want to change the actual "map", as I need many of the widgets that need this tag to load and function correctly.  was thinking of simply loading the mapbox tiles and addinn them to the map.  Following the mapbox API documentation, I loaded the library:

<script src='https://api.tiles.mapbox.com/mapbox.js/v1.6.2/mapbox.js'></script>
<link href='https://api.tiles.mapbox.com/mapbox.js/v1.6.2/mapbox.css' rel='stylesheet' />

I then tried to create the layer of the map box tiles and add that to the map:

var layer = L.mapbox.tileLayer('mapbox.ID');
map.addLayer(layer);

Doing this I get two related errors.

1 - GET http://a.tiles.mapbox.com/v3/mapbox.ID.json 404 (Not Found)
2 - could not load TileJSON at http://a.tiles.mapbox.com/v3/mapbox.ID.json

What am I missing


texts: Merging CSV layers (XYZ files with point coordinates)

texts:
I have list of many XYZ files with points without header an space delimited.
Example:

-877500 -1014000 471.95
-877500 -1013995 472.28
-877500 -1013990 472.61
-877500 -1013985 473.04
-877500 -1013980 473.39

I would like to merge all files into one QGIS point layer. I was trying this tutorial and no luck with Python. Anyway it doesn't describe working with CSV files as I have trouble defining sole VRT with custom delimiter. I couldn't find definitions for them in VRT driver neither CSV driver.
When trying ogrinfo with .xyz file couldn't read the file. If I change the extension to .csv it reads without trouble. I was not aware that <OGRVRTLayer name=" "> has to be the same as CSV file without extension. Finally I managed to combine files with suggestions from comments and using this code:
<OGRVRTDataSource>
    <OGRVRTLayer name="merged">
        <SrcDataSource relativeToVRT=1>merged.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:5514</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
        <Field name="x" type="Integer" src="field_1" />
        <Field name="y" type="Integer" src="field_2" />
        <Field name="z" type="Real" src="field_3" />
    </OGRVRTLayer>
</OGRVRTDataSource>


texts: Stop SQL Server selecting outside buffer

texts:
I am trying to select all the parks within 200m or a railway.  I can easily create and visualise the 200m buffer.  However when i try and find only the parks that intersect it i get all parks returned, not just those in the intersection area.  This is the code i have so far:

SELECT
P.[geogtab] as parks,
P.[id] as park_id
FROM [dbo].[railway]  S
inner JOIN
[dbo].[parks]  P
ON
S.geogtab.STBuffer(200).STIntersects(P.geogtab) = 1

geogtab is my geography column and it is indexed.  Is there an error in my syntax that causes all items from Table P (the parks) ?


texts: Displaying 3D subsoil models in ArcScene

texts:
I was wondering if anybody had any ideas on how to display a 3D subsoil model in ArcScene?

The data we will have will be exported from rockworks in point form and then we need to convert it to TINs. That's what we believe is the correct approach.

I have never done any extensive work in arcscene before and I couldn't find anything too helpful on the internet.


texts: What I can do to find fastly superimpositions?

texts:
I have a project with several geometries belonging to several vectors (more than 10)...
I need to locate the areas where there is a superimposition between two geometries not necessarily belonging to the same vector...

How I can do this operation? I am using QGis 2.12.0


texts: How to calculate polygon area considering slopes?

texts:
I am looking for an accurate algorithm or a service to calculate surface area on earth where points are calculated on the basis of GPS Coordinates.

I am using Google Map Api version 3 and are drawing polygons based on recorded coordinates but I don't think the standard ways of calculating area of the polygon will take into account of slopes(hills). Do I need to work on contours for such thing?

Are there any third party services may be ArcGis or some other that takes into account of slopes as well


texts: .prj for NAD_1983_Canada_Atlas_Lambert not recognized correctly

texts:
I'm a fairly new user of QGIS and I've run into an issue that I can't seem to solve on my own. I'm trying to work with the NRCan Atlas of Canada Waterbodies SHP file but I'm having an issues re-projecting other datasets to line up with this one.

When I reviewed the .PRJ file packaged with the NRCan data, it indicates it is Lambert Conformal Conic projection with two standard parallels at 49°N and 77°N (see .PRJ text below); however, when I add the file to QGIS (v1.8.0-Lisboa) and view the layer properties, the CRS is displayed as follows:

  USER:100001 -  * Generated CRS (+proj=lcc +lat_1=49 +lat_0=49 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs)

Note that the second standard parallel (which should be in there as '+lat_2=77') is missing... I'm hoping this explains why my other datasets aren't lining up when I project them, but I can't seem to get QGIS to recognize the second parallel.

Suggestions? Or is this second parallel not an issue for re-projecting data?

Unmodified NRCan .PRJ file contents

 PROJCS["NAD_1983_Canada_Atlas_Lambert",
 GEOGCS["GCS_North_American_1983",
  DATUM["D_North_American_1983",
  SPHEROID["GRS_1980",6378137.0,298.257222101]],
  PRIMEM["Greenwich",0.0],
  UNIT["Degree",0.0174532925199433]],
 PROJECTION["Lambert_Conformal_Conic"],
  PARAMETER["False_Easting",0.0],
  PARAMETER["False_Northing",0.0],
  PARAMETER["Central_Meridian",-95.0],
  PARAMETER["Standard_Parallel_1",49.0],
  PARAMETER["Standard_Parallel_2",77.0],
  PARAMETER["Latitude_Of_Origin",49.0],
  UNIT["Meter",1.0]]


texts: ArcMap shutting down when I add rasters to a mosaic in foreground processing

texts: I'm trying to add rasters to a mosaic - I was initially trying to do this by background processing but I was receiveing an error - I've followed the advice of this article (http://support.esri.com/en/knowledgebase/techarticles/detail/38695) (turning off background processing/turning off the option to add geoprocessing results to the display) but I still can't get it to work - now instead of returning an error message ArcMap just shuts down with no warning.


texts: Leaflet fire click event by coordinates to specific layer

texts:
I have a county polygon GeoJSON file shown as 4 theme layers, The attributes are the same but the symbology is different. I have a regions file, each region is composed of several counties. I have this as a black outline, so you can see what counties are in what region. So the regions is on top of the turned on themes. What I want to do is get the region to pass the click event to the turned on theme. I'm using radio buttons for the themes so only one can be on at a time. I've tried various variations of this code without luck. The region layer's on-click passes the click X,Y and fires a function that should check the theme layer and use the X,Y for a click.

if (map.hasLayer(theme1)) {
     theme1.fireEvent('click', {latlng:[y,x]});
}
if (map.hasLayer(theme2)) {
    map.theme2.fire('click');
}

Error message is Can't read fire of undefined.


texts:
I ran https://github.com/opengisch/QField build for Android using Powershell with the following steps:

Install docker for Windows
Execute Pull command at https://hub.docker.com/r/opengisch/qfield-sdk/tags/
Run [git clone https://github.com/opengisch/QField.git] to download the source locally
Execute [cd QField] and move to the downloaded folder
Run [git submodule init]
Run [git submodule update]
Run [docker run -v $(pwd):/usr/src/qfield docker.io/opengisch/qfield-sdk:$TAG /usr/src/qfield/scripts/docker-build.sh --rm] , Enter the 2 QField SDK tags in "$TAG"

I get the error

ERROR MISSING APP_VERSION OR APP_VERSION_CODE ENVIRONMENT VARIABLE

What kind of settings are needed to solve this?


texts: Export a table that content results value

texts:
I import a table that have a collection of points defined by (pts), but when I call (pts) by ee.Feature (pts) I can not export the table.

var feat = ee.Feature(pts, {date: ee.Date(date), rain: rainAtDate})

Can I export the table that content (value result (correctFeats) )
for each geometry into the table?

var pts = ee.FeatureCollection('users/mustafalhealy/Z_HdrologyStation')

    // Start and End Dates
    var inidate = ee.Date.fromYMD(2020,01,24)
    var enddate = ee.Date.fromYMD(2020,01,25)

    // Difference between start and end in days
    var difdate = enddate.difference(inidate, 'day')

    // Import GSMaP data
    var gsmap = ee.ImageCollection('JAXA/GPM_L3/GSMaP/v6/operational')
    .filterDate(inidate, enddate)
    .select('hourlyPrecipRate');

    // Time lapse
    var lapse = ee.List.sequence(0, difdate.subtract(1))
    var inidate = ee.Date('2020-01-24')
    var listdates = lapse.map(function(day){
      return inidate.advance(day, 'day')
    })

    // Iterate over the list of dates
    var newft = ee.FeatureCollection(listdates.iterate(function(day, ft) {
      // Cast
      ft = ee.FeatureCollection(ft)
      day = ee.Date(day)

      // Filter the collection in one day (24 images)
      var day_collection = gsmap.filterDate(day, day.advance(1, 'day'))

      // Get the sum of all 24 images into one Image
      var sum = ee.Image(day_collection.sum())

      // map over the ft to set a property over each feature
      var temp_ft = ft.map(function(feat){
        var geom = feat.geometry()
        var value = sum.reduceRegion(ee.Reducer.first(), geom, 1000).get('hourlyPrecipRate')
        return feat.set(day.format('YYYY-MM-dd'), value)
      })

      // Return the FeatureCollection with the new properties set
      return temp_ft
    }, pts))

    // get the dates and the rain intensity per day
    var dates = newft.first().toDictionary().keys();
    var rain = newft.first().toDictionary().values();

    // build a feature collection with property rain and date
    var correctFeats = ee.FeatureCollection(dates.map(function(date){
      var ind = dates.indexOf(date);
      var rainAtDate = rain.get(ind);
      var feat = ee.Feature(pts, {date: ee.Date(date), rain: rainAtDate})
      return feat;
    }));
    print(correctFeats);

    //Making a chart
    //var chart = ui.Chart.feature.byFeature(correctFeats, 'date', 'rain')
    //print(chart);

    // Export the FeatureCollection.
    Export.table.toDrive({
      collection: correctFeats,
      description: 'Rain',
      fileFormat: 'CSV'
    });

Link to demo on Google Earth Engine


texts:
I ceated my table (borne) in PostGIS

id_borne serial NOT NULL,
  num_borne character varying,
  shape_borne geometry

I successfully changed the SRID of my table using

select UpdateGeometrySRID('public', 'borne', 'shape_borne', 26191) ;

The old SRID was 0.

The problem is that when I try to add a new element to my table borne I have this message error:

geometry srid (0) does not match column srid (26191)

That is my problem! The SRID is updated but when I try to add a new point, I get this message.


texts: I have a raster (a catchment area) and a weather station with its coordinates (imported as .txt in ArcGIS) in UTM system. I want to know which row and column of my grid belongs to this station. How can I do it?


texts: Copying attribute values to new layer in QGIS with python

texts:
I am making a buffer for a point layer and want the attribute values to be copied over to the new layer. I am able to duplicate the attribute fields, but not populate them with the values. I get the error
AttributeError: 'QgsVectorDataProvider' object has no attribute 'setAttributes'

lyr = QgsProject.instance().mapLayersByName('lyr')[0]
iface.setActiveLayer(lyr)
lyr = iface.activeLayer()

#set var
width=30

buffLyr = QgsVectorLayer('Polygon?crs=EPSG:3857', ' Buffer', 'memory')
pr= buffLyr.dataProvider()

#duplicate the field from the original layer and update fields
pr.addAttributes(lyr.fields())
buffLyr.updateFields()

#apply buffer to each feature
for selection in lyr.getFeatures():
    buff = selection.geometry().buffer(width_m,6)

    #Add feature attribute values
    attrs =selection.attributes()
    pr.setAttributes(attrs)

    #Create a new feature:
    b = QgsFeature()
    b.setGeometry(buff)

    #Add the feature to the data provider:
    pr.addFeatures([b])

    buffLyr.updateExtents()

QgsProject.instance().addMapLayers([buffLyr])


texts: Controlling how Preset Color Ramps are applied to layers in ArcMap using ArcPy?

texts:
In ArcMap Style Manager, I have a Preset Color Ramp defined. I would like to control how the colors are applied when doing Add All Values, but it seems colors in the set are applied at random and differently in ArcMap and with arcpy.mapping.

When I apply the ramp from ArcMap, Layer Properties, Symbology, Add All Values, these are the results:

When I apply the ramp using arcpy.mapping and lyr.symbology.addAllValues(), these are the results:

It is a totally different set of colors, the orange from the first image and brown from the second are different. I have a Python automated map export process and I would like some consistency in how colors are applied for this layer.

Does anyone have a way to make this more consistent?

Here is the code snippet:

templateMxd = r'D:\...\biodiv_npc_map.mxd'
mxd = arcpy.mapping.MapDocument(templateMxd)
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    lyrList = arcpy.mapping.ListLayers(mxd, "*", df)
    for lyr in lyrList:
        if lyr.name == "Mapped Native Plant Communities near or on property":
            lyr.symbology.addAllValues()


texts:
How can I integrate ArcSDE of Oracle with the QGIS latest version?

I am using ArcGIS 10.4 and QGIS 2.16 version.


texts: mmqgis plugin not showing up in plugin dropdown menu

texts: I have installed the plugin from the plugin dialog box, but I am unable to figure out exactly how to run the plugin.


texts: How to preserve geometry when saving vector files to a SpatiaLite db (QGIS)

texts:
I have been trying to add multiple vector layers to a SpatiaLite database, however using the integrated DB Manager in QGIS 2.14.10 and 2.18.7, whenever I add a shapefile using the "Import Layer/File" button, and follow the menu, the item added does not have geometry, with only the attribute table appearing in the DB.

I can create DBs containing single layers using the "Save as..." right click menu option, but I can't save multiple vector layers into the same DB.

Other users have suggested using the QSpatiaLite plugin, however the plugin doesn't run when installed, and hasn't been updated in 6 years. QspatiaLite seemed like a good answer, if anyone knows of a similar plugin that has been upkept, please let me know.

Here is a snapshot of my SpatiaLite DB, with GLA_Zips being the vector file that was added to the original combined_mix sqlite file.  Despite checking the geometry column option in the import menu, the geometry still isn't imported.  When I click add to canvas in the context menu for the item, its added to the layers panel as a table, not a polygon vector file.

The original combined_mix layer looks like this in the DB, for comparison, and works in the layers panel just fine.


texts:
This is the error I am getting in QGIS Bonn 3.2.2:

Traceback (most recent call last):
  File "/usr/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 3, in <module>
NameError: name 'QgsMapLayerRegistry' is not defined

This is the script i am trying to:

from qgis.PyQt.QtCore import QVariant

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
    if layer.wkbType() == QGis.WKBPolygon:
        with edit(layer):
            layer.dataProvider().addAttributes( [ QgsField("Node_ID", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("Node_Type", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("T_1", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("R_1", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("S_1", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("QS_1", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("T_2", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("R_2", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("S_2", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("QS_2", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("T_3", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("R_3", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("S_3", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("QR_3", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("T_4", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("R_4", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("S_3", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("QR_4", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("Meridian", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("Disc_Mon", QVariant.String) ] )
            layer.dataProvider().addAttributes( [ QgsField("Dis_Lat", QVariant.Double) ] )
            layer.dataProvider().addAttributes( [ QgsField("Dis_Long", QVariant.Double) ] )
            layer.dataProvider().addAttributes( [ QgsField("Dis_Dist", QVariant.Int) ] )
            layer.dataProvider().addAttributes( [ QgsField("Dis_Azi", QVariant.Int) ] )
    else:
        pass


texts: Comprehensive database of all cities/towns in South Africa?

texts:
The title says it all really - I'm looking for a comprehensive list/database of all cities/towns in South Africa. I know about these guys: http://www.afrigis.co.za/, but they're just not getting back to me no matter how many times I contact them. They were recommended to me as the only really trustworthy and comprehensive source for the data I need, but they're obviously not extremely helpful.

If anyone has any info on where I can get this then please let me know. Yes, I have tried Google and no I haven't found anything that is comprehensive enough (except for AfriGIS).


texts:
I have Trimble Juno 3 GPS and I can't go into TerraSync.

The GPS does not receive signals from satellites.

This happened suddenly after I deleted all the files in the file manager and then I turned it off.


texts:
I have multiple shapefiles in a folder which I am trying to compare with same named features in geodatabase. I want to first clip within an area polygon and then copy and overwrite the features in the geodatabase. So far the describe for shapefile names isn't working;

# Import arcpy module
import arcpy

# Set workspace

arcpy.env.overwriteOutput = True
shpFile = arcpy.env.workspace
arcpy.AddMessage("\n")

shpName = arcpy.Describe(shpFile).baseName
featureName = arcpy.Describe(geoPath).baseName

# Count the number of shapefiles clipped
count1 = 0
# Count the number of feature classes being copied
count2 = 0

try:
 # Loop through shapefiles and clip to region
 for shp in shpFile:
     if shpName == featureName:
         arcpy.AddMessage(shpFile)
         arcpy.Clip_analysis (shpFile, clipRegion, outFc, "", "", "")
         count1 += 1
         # Copy clipped shapefile into file geodatabase
         arcpy.FeatureClassToFeatureClass_conversion(outFc, geoPath + "\\" + shpName)
         count1 += 1
     else:
         arcpy.AddMessage("Error: Shapefile/s not clipped and copied to file geodatabase ")


texts: Can I create extension for mapd database?

texts:
I want to create extension for mapd database to access from geoserver.

so I want to make sure mapd support GeoServer or not?


texts: Numeric scalebars don't (always) update when changing the map's extent in QGIS

texts:
If you change your map's extent but maintain the scale in a QGIS Layout, the scalebar's width will be updated proportionally with the CRS's distortion for the new extent. E.g. if your map (crs 3857) is focussing on the equator, and you drag the map's extent to the north (maintaining the same scale), the scalebar's width will enlarge automatically, which is good, because in pseudomercator areas around the equator are displayed relatively smaller than areas futher away from the equator.
However, if you do the same thing while having a 'numeric' scalebar, the scale is not updated in the same way. Why is that?

And how come, if you do the same thing while your map has a geographic CRS such as WGS84 (4326), the numeric scale will be updated just like the scalebars are?


texts: How do I restrict a fme transformer to work only within a polygon mask?

texts: I would like to define the area in which a transformer should be applied by using a polygon mask. I know this can be done in FME but I can´t find any references to it. Can anyone enlighten me?


texts: I am getting error in PostGIS 'function st_clip(raster, integer, geometry) does not exist'?

texts:
I have to find the count of the number of pixels in a given band of a clipped raster.
For clipping, I am using ST_clip(). Through clip function raster is clipped with a polygon.
 I am getting following error

  function st_clip(raster, integer, geometry) does not exist

I am not sure about the placements of the paranthesis.

SELECT ST_ValueCount(ST_Clip(rast,1,ST_Transform(ST_GeomFromText('POLYGON((76.36792879303643 13.633268921278498,77.64804272121083 13.625811445922459,77.58727436708882 12.839858067425443,76.54887536367008 12.766426303994999,76.36792879303643 13.633268921278498))',4326), 4326)),1,false) AS valuecount FROM raster_file;

EDITED

Should I use ST_UNION before ST_clip as my raster file is tiled and then saved in postgis using 'raster2pgsql' command.

Otherwise please provide a proper example to clip a raster from polygon, and to find the valuecount.


texts:
I want arcpy to perform the following: load a shapefile, prepare to edit it, open its attribute table, perform a numpy.digitize using the values of a certain attribute and some bins provided separately.
My  code looks like this:

network = r"C:\Mylocation\MyShapefile.shp"

intersected = arcpy.da.UpdateCursor(network, ["DIAMETER", "Other_Attributes_I_Will_Need_Later"])

bins = [ bin1, bin2, bin3, ... binn]

cats = numpy.digitize(intersected.DIAMETER, bins)

But this way I receive the error "AttributeError: 'str' object has no attribute 'DIAMETER'. I guess the problem is that it my code does not "dig" into intersected to find the values I need, but therefore how can I solve this problem?


texts: How can I quantify whether locations of one 'type' are closer to one another compared to a second 'type'?

texts:
This might have an obvious answer but I'm new to GIS and can't seem to find the right technique... I've got GPS points for a bunch of households and have calculated the distance (in km) between each pair of houses. Each household can be one of two types (which I have assigned) - let's say 'jock' and 'nerd'. I want to quantify whether nerd households are physically closer to one another within a neighbourhood compared to jock households. I'm kind of looking for a measure of isolation or segregation... What would be the best approach?

(NB: there are not equal numbers of jock and nerd households in the dataset.)

I've looked into techniques such as cluster analysis, principle coordinate analysis and dissimilarity matrices, but none of these seem right. I've even tried doing t-tests and Mann-Whitney tests on the distances categorised by household type but there are so many dyads that the tiniest difference between groups comes out as statistically significant.


texts: Exporting selected features using ArcPy?

texts:
# Wetlands buffer.py
# Description: Program that determines if Wetlands are within 0.5 miles of specific point files with the program returning a 'yes' or 'no' if a wetland is within the range or not.

# Import ArcGIS python library
import arcpy

# Define workspace
arcpy.env.workspace = "P:\Python"

# create buffer around single substation
substation_buffer = arcpy.Buffer_analysis("substation.shp", "P:\Python\substation_buffertest", "0.5 Miles", "FULL", "ROUND", "LIST", "", "PLANAR")
print "Made buffer layer"

# turn substation_buffertest into a feature layer for the select by location tool
fl = arcpy.MakeFeatureLayer_management("substation_buffertest.shp", "buffer_layer")
print "Made feature layer"

# select wetlands within buffer
wetland_selection = arcpy.SelectLayerByLocation_management("buffer_layer", "CONTAINS", "Wetlands.shp", "", "NEW_SELECTION", "NOT_INVERT")
print "Selected by location"

I'm trying to write a program for the company I work for that determines if a point file's buffer zone contains wetlands from another feature class. I'm trying to do this through a search cursor but the attribute table that is created only contains the FID, Shape, and Id. I've made it this far but I'm not so sure where to go from here.


texts: Associating location data from PostGIS for a Choropleth map?

texts:
I'm trying to build my first choropleth map from data that we have stored in our PostGIS enabled Postgres database.

So far I've been able to build a map on Mapbox using a topojson set of UK postal area data, but I'm stumped on how to show concentration of our data points per postal area from our database.

Basically I'm trying to show we have a concentration of user type X in postal areas.

I know enough that querying the db point-in-polygon for each polygon in the postal data set would be slow and cumbersome, but is there a way to search the data by postal area and return a merged data set e.g. x users in area 1, y users in area 2?

Should I import the postal areas into a spatial table to do the query against in PostGIS?

Currently all the data points are output as JSON but I'm not sure how to associated the relevant points with each postal area.


texts: Error when building network dataset in ArcGIS Pro?

texts:
I have a network dataset with includes roads and turn feature class.
I updated it once a week. I use ArcGIS Pro to create the network dataset based on a template (but can do manually with network wizard and all errors are the same).

I import turn feature class from the previous network dataset and roads comes from somewhere else.

After build is done, I always get this error for all turns:

  No edge element can be found corresponding to the turn identifier

I don't understand what this error means and how to fix it.
If I move all turns manually slightly, then build, it works. But it won't work for me as I need automate the process of road network update.


texts: ESRI 10.5 Azure Cloud builder install - LoadBalancer

texts:
After installing ArcGIS Enterprise (aka ArcGIS for Server 10.5) an Azure LoadBalancer instance is created, in addition to a configured ESRI webadaptor being installed.  It's unclear (not documented) what role the LoadBalancer plays and how it works with the ESRI web adaptor.  It seems that they potentially do the same thing which is load balancing.  Below is a screen shot of a single server install.  In a multi-site install the contents of the resource group are the same with the exception that there are more hosts.  The creation of the LoadBalancer is new to 10.5, what is it doing?

This question is covered for AWS and can be found on an ESRI FAQ titled Is ArcGIS Web Adaptor installed when I create my site with Cloud Builder? but a similar answer for Azure does not exist.


texts: Shooting star algorithm in postgres 8.4 and postgis1.4 showing fatal error and remaining are working fine

texts:
When I try to run this shooting star algorithm, the postgis suddenly loosing connection to postgres and showing error no connection to server

SELECT * FROM shortest_path_shooting_star('
   SELECT gid AS id,
      start_id::int4 AS source,
      end_id::int4 AS target,
      shape_leng::float8 AS cost,
      x1, y1, x2, y2,
      rule, to_cost
   FROM viz_network',
   173,
   97,
   false,
   false);

And digkstra algorithm is working fine :

SELECT * FROM shortest_path('
  SELECT gid AS id,
      start_id::int4 AS source,
      end_id::int4 AS target,
      shape_leng::float8 AS cost
  FROM viz_network',
  1,
  310,
  false,
  false);


texts:
I'd like to save the R-tree spatial index generated by GeoPandas to disk.
GeoPandas automatically generates an index for each shape (points in this case) in the dataframe.
I'm attempting to create a new index using the R-trees built in built-in disk serialization, but not sure how to build it up to match the original dataframe.
I've made a hack by simply iterating over the number of shapes and adding bounds... this seems to work for my purposes but am quite sure this is not the correct approach.
Apologies for not taking the time to fully understanding how R-trees work.
I borrowed from this post Saving python Rtree spatial index to file? and found this discussion https://github.com/geopandas/geopandas/issues/426
import numpy as np, pandas as pd, geopandas as gpd
from shapely.geometry import Point
from rtree.index import Index

lats = np.arange(10, 20, 1)
lons = np.arange(60, 80, 2)
df        = pd.DataFrame({'lat': lats, 'lon': lons})
df['geo'] = df.apply(lambda df: Point(df['lon'], df['lat']), axis=1)
gdf       = gpd.GeoDataFrame(df, geometry='geo', crs='EPSG:4326')
sidx      = gdf.sindex
bounds    = sidx.bounds

idx       = Index('path')
for sid in range(sidx.get_size()):
    idx.insert(1, sidx.bounds)
idx.close()


texts: Writing ArcPy Update Cursor "Where Clause" that Refers to its own Rows?

texts:
I'm using Arc 10.2.2.

I have 2 files:

"Input_Table.exe" (only 3 columns: "ID", "Lat", and "Lon")
"Target_Shapefile.shp"  (many columns including the same "ID", "Lat", and "Lon")

Any rows in the Target_Shapefile with an ID that's present in the Input_Table need to be updated with those appropriate values.

I'm aware of the Update Cursor and that using its 3rd parameter, the "Where Clause", will reduce the number of retrieved rows to be updated, thus saving time.

But how do I refer to a row in the same shapefile that's being updated?  (Target_Shapefile)

More specifically, is the following syntax for my "Where Clause" valid?

    Where = "Where ['ID'] in Input_ID_list"

My Complete Code:

# Create a list of all IDs to be updated in the Target_Shapefile.
rows = arcpy.da.SearchCursor(Input_Table)
for row in rows:
    Input_ID = row.getValue('ID')
    Input_ID_list = Input_ID_list + Input_ID

# Set "Where Clause"
Where = "Where ['ID'] in Input_ID_list"

# Create the Update Cursor...
with arcpy.da.UpdateCursor(Target_Shapefile,['ID', 'Lat', 'Lon'] , Where) as Tcursor:
        for row in Tcursor:
                do some things...


texts: Do you know about any site where there is any discussion on Geomedia products (besides Intergraph Synergy site which is not very active lately)?


texts:
I'm new to PostGIS and a bit lost...
I work with linestrings that sometimes overlap, for example:

I'd like to distinguish these linestrings by splitting them like Google Maps do:

I was thinking about using ST_Intersection() and ST_OffsetCurve(), but maybe the problem is more complex than i thought ?
Or maybe there is a simpler solution i did not see ?


texts: Geopandas simplify results in gaps between polygons

texts:
The goal is to simplify the geometries of a shapefile containing multipolygons. However, when I use Shapely's simplify in geopandas, the result contains gaps between the polygons. I was hoping that preserve_topology would avoid this? How can I avoid or fix the issue?

gdf_simplified = gdf.copy()
gdf_simplified["geometry"] = gdf.geometry.simplify(tolerance=TOLERANCE,preserve_topology=True)

my geoDataFrame is in in {'init': 'epsg:4326'}and I use a tolerance of 360/43200 which corresponds to 30" degrees (appr. 1 km at the equator).

The old polygons are displayed using the black line. The new resulting polygons are depicted using colors.

For now, I was able to revert to MapShaper and get simplified shapes without the gaps. However, I would prefer a python based approach.


texts: mapcache_seed - "No space left on device"

texts:
Introduction

I'm trying to seed a cache for a service, that already exists using a second grid.

The error seems to be a problem with the file system, not the seed itself. If I do not create a seed, but request a tile manually I get the exact same error.

Information on the operating system

MapCache and the mapcache_seed command are executed on a Debian Stretch Box
The INTERNAL_ADDRESS refers to another local Debian Stretch Box with an out of the Box MapServer installed
The actual Data of the Map comes from another local Box, that runs PostgreSQL 9.4 / PostGIS on Debian Jessie

MapCache is running on a Debian Stretch Box

Information on the file system

The system has 60 GB of space mounted on /mapcachedata. The space needed for one cached service is about 50 GB. To add another cache I added another 60 GB on /mapcachedata/cache/karte_sl/UTM32. I believe the problem might have to do with mapcache not handling this correctly. I inherited the system and tried to add this cache with minimal effort. The old files are stored in /mapcachedata/cache/karte_sl/GK2. Creating a few GB of files in /mapcachedata/cache/karte_sl/UTM32 effected the correct hard drive.

Error Messages

The command used was mapcache_seed -c /opt/mapcache/mapcache.xml -t karte_sl -g UTM32. The console response was:

failed to seed tile z0,x0,y0:
failed to create directory /mapcachedata/cache/karte_sl/UTM32/00/000/000/000/000/000/000.jpg: No space left on device

aborting seed as 100.0% of the last 1000 requests failed

failed to seed tile z1,x0,y0:
failed to create directory /mapcachedata/cache/karte_sl/UTM32/01/000/000/000/000/000/000.jpg: No space left on device

aborting seed as 100.0% of the last 1000 requests failed

The XML File looks like this: https://www.codepile.net/pile/dZ3z51qP

When I request a tile on via HTTP GET /mapcache/tms/1.0.0/karte_sl@UTM32/4/2/0.jpeg I get pretty much the same response:

failed to create directory /mapcachedata/cache/karte_sl/UTM32/04/000/000/000/000/000/000.jpg: No space left on device

Question

Does it have to do with the mounted hard drive, or is there another problem?


texts:
i'm trying to import raster to postgis using this command ;

raster2pgsql -I -C -s 4326 -t 100x100 -M "C:\DATA TRIAL APOLLO\mosaik129_130.ecw" -F public.mosaik129_130 | psql -d apollo_db

but then error message :

ERROR: unable to read raster file: C:\DATA TRIAL APOLLO\mosaik129_130.ecw

and it's requiring a password, even on my windows account, i'm not setting any password.

does anyone have solutions?


texts: Converting polygon to raster using ArcGIS Desktop?

texts:
I'm just learning the basics of Arcmap and I'm currently stuck on doing this.

Ive got a road with a 500m buffer (polygon) and I'm attempting to convert it into a raster so I can reclassify it into (1(suitable) and 0(unsuitable) in order to use the data for the boolean tool so I can have 2 factors merge into 1 area.

However I'm currently stuck on converting the polygon into actual data that I can reclassify because currently I'm just getting an errors, or no data.


texts:
I'm trying to build 3D buildings from 2D footprints. I extrude by the building heights, and input the surface by editing the properties. From here I use the Layer 3D to Feature Class tool to convert to multipatch.

However the output from this results in a seemingly empty feature class, with no entries in the Attribute Table. I've tried saving the extruded buildings as a layer and using that as input to the Layer 3D to Feature Class tool, producing the same result of an empty feature class.

Am I missing something about this tool? Thanks for any help.


texts: Theory notions on exporting bufr data in geotiff format

texts:
I'm trying to export some bufr data from polar satellites/sensors like DMSP, ASCAT, ATMS etc in geotiff image.
When my program rasterizes satellite data (brigthness temperature for example) by converting long/lat coords in pixel coords, I obtain an image with scattered data. I think it's because satellite data are sampled in a not equispaced manner due to the Earth curvature, and then I map this data to an equispaced matrix of points (my tiff).
My question is simple: is my theory right?
I think yes but some online readings put me in doubt.


texts: How to resolve pgrouting dijkstra_sp() function error?

texts:
insert into dijsktra_results(geom) select ST_UNION(the_geom) from dijkstra_sp('roads_layer',
            (select source_id from roads_layer where gid=9975)::integer,
            (select target_id from roads_layer where gid=2967)::integer);

This query to perform the routing operations.
"dijsktra_results" is a newely created table for routing operations with gid and geom columns of 2D(MultiLineString).
This is the error I am getting

ERROR:  column "the_geom" does not exist
LINE 1: SELECT gid,the_geom FROM shortest_path('SELECT gid as id, so...
                   ^
QUERY:  SELECT gid,the_geom FROM shortest_path('SELECT gid as id, source::integer, target::integer, length::double precision as cost FROM roads_layer', '1' , '20142' , false, false), roads_layer where edge_id = gid
CONTEXT:  PL/pgSQL function "dijkstra_sp" line 13 at FOR over EXECUTE statement

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

ERROR: column "the_geom" does not exist
SQL state: 42703
Context: PL/pgSQL function "dijkstra_sp" line 13 at FOR over EXECUTE statement

Can anyone resolve this problem or guide for an extent am I gone wrong anywhere.Till now everything worked fine.I calculated vertex_id,edge_id,length....
And how this, dijkstra_sp() function works and calculates the shortest path.


texts:
I'm trying to setup a process using PyWPS. I've been able to send geometry and receive them in GML with this code :

     self.vectorin = self.addComplexInput(identifier="vectorin",
                                              title="Geometrie 1",
                                              formats=[{"mimeType":"text/xml"},{"mimeType":"application/xml"}])
     self.vectorin2 = self.addComplexInput(identifier="vectorin2",
                                              title="Geometrie 2",
                                              formats=[{"mimeType":"text/xml"},{"mimeType":"application/xml"}])
     self.vectorout = self.addComplexOutput(identifier="vectorout",
                                              title="Vector file",
                                              formats = [{"mimeType":"text/xml"}])

 def execute(self):
    self.vectorout.setValue(self.vectorin.getValue())
    return

This operation is useless but it is the only exemple with PyWPS 3.2.x using ComplexData. I'm trying to extract the geometry from both of my inputs to execute operations with GDAL/OGR.

Is there anyway I can do it and then send it back in GML as the output?

Update :

Exemple : when I'm trying to run

def execute(self):
  geometry1 = Osgeo.ogr.CreateGeometryFromGML(self.vectorin.getValue()).GetGeometryRef()
  geometry2 = osgeo.ogr.CreateGeometryFromGML(self.vectorin2.getValue()).GetGeometryRef()
  geometry3 = geometry1.Intersection(geometry2)
  geometry3.SetGeometryDirectly(geometry3)
  self.vectorout.setValue(geometry3.ExportToGML())
  return

I get "Failed to execute WPS process [testPerfoPy]: 'NoneType' object has no attribute 'GetGeometryRef'" So I'm presuming that I can't use CreateGeometryFromGML ? I can't stop the code while it executes so I'm not even sure if the value in self.vectorin.getValue() is really the GML.


texts: Efficiency in PostGIS: objects near point- ST_Intersects around buffer Vs ST_DWithin

texts:
I need to find all buildings in radius of 25m from given location. I thought about 2 PostGis quires that can handle this, when I run it they both give same time execution. But I want to know which one of them is the more efficient way?

SELECT *
FROM buildings
WHERE
ST_Intersects(ST_Buffer(ST_Transform(ST_GeomFromText('POINT(lon lat)', 4326), 2163), 25), ST_Transform(buildings.geom, 2163));
-----------
SELECT *
FROM buildings
WHERE
ST_DWithin(ST_Transform(ST_GeomFromText('POINT(lon lat)', 4326), 2163), ST_Transform(buildings.geom, 2163), 25);


texts: Separate boundary from adjacent polygons

texts:
I have a list of polygon features that share a common boundary. Is it possible to extract that line somehow? I am thinking about making a spatial join that shows me polygons that share a line segment, but how can I extract the line that is the reason for their join.
For example, in the picture, polygons 454 and 218, 1147 and 684, 684 and 56 share a boundary (circled in blue).


texts: Record distance moved and new lat/long in attribute table after snapping points to line

texts: I am using QGIS 3.6. I have a dataset of fish locations (points) that I snapped to a river (line). I would like to add the distance each point moved to the attributes table, as well as the new lat/long. As there are several hundred fish localities, I would like to do this as a batch after snapping them all rather than manually measure the distances between each pair of points.


texts:
Is there a package to call QGIS functions directly from R?
Is rPython the best alternative?


texts: Using QGIS to crop shapefiles using raster image

texts:
I have a topographic map in raster file and I want to add the lines for fault systems on top (in shapefile), but I want to crop them so I'm only left with those that are within the raster image.

Is there a way to use the raster to select the shapefile that are within the boundaries of the raster file?

I'm using QGIS 2.10.1.


texts: Displaying tiles build with gdal2tiles in QGIS 2.14

texts:
I'm trying to create a WMTS service showing an orthophotography.
I've build the tiles with gdal2tiles (https://gdal.gloobe.org/gdal/gdal2tiles.html) and it seems to be ok as i can see the tile in the viewer (leaflet.html) automatically generate by

gdal2tiles -s EPSG:2154 -z 16-21 Z:/01_A43/04_RASTER/ORTHO/952500-6500500
.tif C:/ORTHO_TILES

The result of such command is a repository (ORTHO_TILES) containing sub-directories (Zoom Level, 16 --> 21) containing sub--directories (referred to X in the world grid), containing tiles which have the name of the Y reference in the world grid.

example : C:/ORTHO_TILES/18/135615/168420.png

looks like this

then i try to visualise this in QGIS with an .xml file. I wrote it as it's specified in http://www.gdal.org/frmt_wms.html

 <GDAL_WMS>
  <Service name="TMS">
    <ServerUrl>file://C:/ORTHO_TILES/${z}/${x}/${y}.png</ServerUrl>
  </Service>
  <DataWindow>
   <UpperLeftX>-20037508.34</UpperLeftX>
   <UpperLeftY>20037508.34</UpperLeftY>
   <LowerRightX>20037508.34</LowerRightX>
   <LowerRightY>-20037508.34</LowerRightY>
   <TileLevel>21</TileLevel>
   <TileCountX>1</TileCountX>
   <TileCountY>1</TileCountY>
   <YOrigin>top</YOrigin>
  </DataWindow>
  <Projection>EPSG:3857</Projection>
  <BlockSizeX>256</BlockSizeX>
  <BlockSizeY>256</BlockSizeY>
  <BandsCount>4</BandsCount>
  <Cache />
</GDAL_WMS>

The layer is load in QGIS but nothing appears in the canvas...

In the past, i'va done similary tricks using QTiles to create the TMS repositories and it worked well. I don't understand why with gdal2tiles it's different.

As QGIS doesn't push any error message and a I can zoom on the extent of the layer, it seems that QGIS recognize the xml as a raster without issues.

i have already read this question but there is no answer that worked for me here

gdal2tiles and qgis


texts:
I am using leaflet map. I have copied the map text from one of it's example.

var mymap = L.map('mapid').setView([51.505, -0.09], 13);

L.tileLayer('https://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/{z}/{y}/{x}', {
    attribution: 'Tiles &copy; Esri &mdash; Esri, DeLorme, NAVTEQ, TomTom, Intermap, iPC, USGS, FAO, NPS, NRCAN, GeoBase, Kadaster NL, Ordnance Survey, Esri Japan, METI, Esri China (Hong Kong), and the GIS User Community',
    maxZoom: 18
}).addTo(mymap);

Is it legal to reduce copyright text ?

Current copyright text

Tiles © Esri — Esri, DeLorme, NAVTEQ, TomTom, Intermap, iPC, USGS, FAO, NPS, NRCAN, GeoBase, Kadaster NL, Ordnance Survey, Esri Japan, METI, Esri China (Hong Kong), and the GIS User Community

I want to change it to

© OpenStreetMap contributors

Edit

I have udpated the map tile provider to openstreet maps.

Reference - http://server.arcgisonline.com/arcgis/rest/services


texts: Writing less than or equal in label of graduated classes in symbology of ArcMap?

texts:
I'm using ArcGIS 10.3.

How do I write less than or equal in the break values of the labeling in symbology?


texts: Generate multiple contour lines at specific elevations

texts: Can anyone suggest a method by which I can produce several user-defined contour lines (from a DEM) within a single shapefile using QGIS/GRASS (2.12.x) and its plugins. Each approach I have tried has only allowed my to create contours in specific intervals, rather than at specific elevations


texts: Is it possible to join two different data sets in QGIS?

texts:
I have one data set which contains all federal states of Germany. The other dataset is about the population density in every km² in the EU.

My goal is now to join these two Datasets so I can answer different questions like 'how many km² with more than 5000 people are in federal state x?'

The first data set looks like this:

and the second like this:

.

I am very new to GIS programs, but from my intuition it sounds like a very basic need to combine different datasets in order to answer a question.


texts: Geo-referenced raster in QGIS does not align with shapefile

texts:
I am using georeferencer in QGIS to align a rater with a boundary shapefile, but after georeferencing, the resulting raster shifts a bit but does not shift completely to the GCP as it is supported to be. As illustrated in the image below, the circled area is supposed to be stretched to align with the red line above it but now there is a gap there.
This problem goes away when using ArcGIS, but there is another problem with exported raster from ArcGIS so I prefer to do it in QGIS.
Does anyone know how to fix it?
georeferenced raster in QGIS:

georeferenced raster in ArcGIS:


texts: How to extrude polygon features from a TIN surface in ArcScene?

texts: How do I Extrude polygon feature class from a TIN surface using ArcScene 10.1?


texts:
I have OS street view mapping as raster images. Is it possible (without the obvious of tracing each dwelling) to identify all polygons on each tile.

To clarify,

My ultimate aim is to produce a shapefile layer of polygons for each component....particularly dwellings in order to create a 260m buffer around each which is easy done.

Using QGIS


texts: Visualizing information when click in polygon - Spatial join

texts:
I am new to cartodb, right now I am trying to do something similar to the cartodb tutorial: Create your first map with CartoDB.js, but with polygons using my own map.

For example:

I did a spatial join with two tables, one has the coordinates and the other one the polygons, once I do the spatial join, three or two polygons appear in the screen of the map, one of them containing 10 rows from the table and the other 20 rows. I want to visualize the count based on this after I click the polygon.

I was wondering if you could help me with some information or tutorial to be able to visualize it?

Thanks a lot.

Best regards,

Walter


texts: List index out of bounds [-1] when making bioclimatic variables in DIVA GIS

texts:
I am working with climate data and I would like to make 19 bioclimatic variables which I can use for ecolgical niche modeling in Maxent. I am following this tutorial http://www.ccafs-climate.org/downloads/docs/Ramirez_Bueno-Cabrera_2009_tutorial_bcvars_creation.pdf (google: bioclimatic variables + Maxent)

After I create the CLM files in DIVA-GIS, I go to data>climate>map and I receive this error:
"List index out of bounds [-1]".  Has anyone run into such an error and do you know the solution?