texts: Intersect does not remove features, leaves them empty. Python

texts:
So I have a shapefile with 16k linestrings, and one polygon of a smaller area.

My intention is to do intersection and keep only the lines or fragments of lines that lie inside the polygon, following this tutorial: http://geopandas.org/geometric_manipulations.html

This is the output of the geometry of the large shapefile:

> network
0        LINESTRING Z (524669.5552664882 6066308.874914
...
16047    LINESTRING Z (735457.2357274202 5693634.842955...
Name: geometry, Length: 16048, dtype: object

this the polygon:

> shtest_utm geometry 0  POLYGON ((781557.9314549405 5812460.808320792,...

This is my code:

file = 'test.geojson'
with open(file) as f:
    data = geojson.load(f)
vecfilename = 'strecken_polyline.shp'
pdb.set_trace()
vector = GeoDataFrame.from_file(vecfilename)
vector.crs = {'init': 'epsg:4326', 'no_defs': True}
vector_utm = vector.to_crs(epsg = 32632)
shtest = GeoDataFrame.from_features(data["features"])
shtest.crs = {'init': 'epsg:4326', 'no_defs': True}
shtest_utm = shtest.to_crs(epsg = 32632)
network = vector_utm['geometry']
network.buffer(20)
network.plot()
shtest['geometry'].plot()
inter = network.intersection(shtest_utm)

And the output of inter is:

  inter
  0        ()
  ...
  16047    ()
  Length: 16048, dtype: object

Is this a normal behavior? I would not expect so. Like that I can not see if there is any line remaining.

Furthermore, .plot() does not work, but that is another problem.


texts: How to define CRS (WGS84) of a mosaic raster file while all the parental raster files are in WGS84 in Python?

texts:
out_meta.update({"driver": "GTiff","height": mosaic.shape[1],
             "width": mosaic.shape[2],
             "transform": out_trans,
             "crs": "+ "+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs"})

Using this code lead me wrong to define the mosaic CRS as the mosaic raster suppose to be in WGS84. After defining the mosaic in WGS84 I will project the file into a local coordinate system. How can I define a custom CRS?
I am totally a new beginner in Python. I am seeking an answer with examples.


texts: Display encoded polyline with OpenLayers

texts:
I try to display routing path on OpenLayers map, but the path is rendering wrong. Back-end returned encoded polylines are displayed in Africa, but it should be in Poland (image below).

I tried to modify example code from OpenLayers map website (http://openlayers.org/en/latest/examples/feature-move-animation.html).

My code:

import Feature from 'ol/Feature.js';
import Map from 'ol/Map';
import View from 'ol/View';
import XYZ from 'ol/source/XYZ';
import Polyline from 'ol/format/Polyline';
import Point from 'ol/geom/Point.js';
import VectorSource from 'ol/source/Vector.js';
import {Group as LayerGroup, Tile as TileLayer, Vector as VectorLayer} from 'ol/layer.js';
import {Circle as CircleStyle, Fill, Icon, Stroke, Style} from 'ol/style.js';
import Projection from 'ol/proj/Projection';
      var polyline = [
            "wcauHspjbBBC@G@G?GAGAECECAG?I@m@HE@",
        "cfauHcrjbBuPoEc\\_GiSiOhCg^bNea@pb@bK|b@hRx^|Ndj@bBzO~Sz]ff@`YhRpYqG",
        "y~ztHayibBKIIIa@]EAC@]x@ABGPMVGPABUl@M^Sj@Qv@MIOIGCECVwAFe@Di@DyA@]@k@?C",
        "{d{tHo|ibB}F~iAg\\`d@_]jV}b@eI",
        "_l~tHo~ebBl@?n@A`Ca@bAAt@l@bAl@Ct@R|KBjAFzAAhB@jA"
         ].join('');

        var route = /** @type {module:ol/geom/LineString~LineString} */ (new Polyline({
            factor: 1e6
        }).readGeometry(polyline, {
            dataProjection: 'EPSG:4326',
            featureProjection: 'EPSG:3857'
        }));

        console.log(route);

        var routeCoords = route.getCoordinates();
        var routeLength = routeCoords.length;

        var routeFeature = new Feature({
            type: 'route',
            geometry: route
        });

        var styles = {
            'route': new Style({
            stroke: new Stroke({
                width: 6, color: [24, 170, 13, 0.8]
            })
            }),
            'icon': new Style({
                image: new Icon({
                    anchor: [0.5, 1],
                    src: 'data/icon.png'
                })
            }),
            'geoMarker': new Style({
                image: new CircleStyle({
                    radius: 7,
                    fill: new Fill({color: 'black'}),
                    stroke: new Stroke({
                    color: 'white', width: 2
                    })
                })
            })
        };

        var vectorLayer = new VectorLayer({
            source: new VectorSource({
                features: [routeFeature]
            }),
            style: function(feature) {
                return styles[feature.get('type')];
            }
        });

        this.map.addLayer(vectorLayer);

        var geoMarker = new Feature({
            type: 'geoMarker',
            geometry: new Point(routeCoords[0])
        });


texts: Displaying multiple images from same point in QGIS

texts:
After the dead of eVis plugin I have troubles working on my project in QGIS 3.16:
I have a shapefile with multiple row, each one with the same position: the content is data from archaeological sites and is derived from a query on a spatialite database.
Each row is provided with the path of a picture.
I need to display, in sequence, all the pictures related to the site (with also their data, stored in the single photo row, but that is ok). Until now the eVis browser from eVis plugin was doing exactly this, but now is dead and I cannot find any replacement.
I've tried to work on the Attribute Forms in the layer properties, but I can only manage to do the required task (via Identify Features) when there is only one image related to a single position.
With multiple images the Identify Features open a tree form, with the path of the images but not the pictures. I attach two example screenshots from a test project (single and multi image point).
Any idea on how to resolve the problem?


texts: How to load OpenStreatMap dump to PostgreSQL?

texts:
I found tool named imposm, I want to make local copy of OSM data and load them to PostgreSQL.

I can't understand how to do it. I run it with follow command:
imposm3 run -config /apps/imposm/config.json

Config is next:

{
    "cachedir": "/home/osmdata/diff/",
    "connection": "postgis://postgres:123456@localhost:5432/geodb?sslmode=disable&prefix=NONE",
    "mapping": "/apps/imposm/mymap.json",
    "srid": 4326
}

I am getting error:
[diff] unable to read last.state.txtstat /home/osmdata/diff/last.state.txt: no such file or directory

What should include this file. If I do not want any "last state" and doing the fresh import.


texts:
In ogr2ogr or ogrinfo, I can successfully query an Oracle database using this command as a test:

ogrinfo -al -so "OCI:USER/PASSWORD@SERVER/DATABASE" -sql "SELECT * FROM SCHEMA.TABLE"

However, building up the query and reusing older queries, I need to get the -sql flag input from a text file:

ogrinfo -al -so "OCI:USER/PASSWORD@SERVER/DATABASE:SCHEMA.TABLE" -sql C:/TEMP/sql.sql

The return this answer from the database: ERROR 1: ORA-00900
The sql.sql file holds same query (SELECT * FROM SCHEMA.TABLE) as above. I've tried to change the sql.sql files quotes with/without/doubles and others with no luck.

How can I format the sql text file, or ogr flag, so it would be read correctly by Oracle 11g.


texts: Use vector layer attributes as edge properties of a QgsGraphArc with PyQGIS

texts:
With PyQGIS 2.18 I created a QgsGraph as describe in the PyQGIS developer cookbook, Release 2.18, ch. 21.2, pp. 95-97. In the cookbook it says

"Vector layer attributes and length of an edge can be used as the properties of an edge".

By now, my source vector layer has variety of attributes (road names, landuse, ...) I would like to use as edge properties for styling purposes, labeling etc of my network analysis results.

I tried to access the edge properties by QgsGraphArc.properties() but this gives my a one-element list of the cost attribute used for graph analysis.

So the Question is: How to use the attributes of the source vector layer for a QgsGraph as properties of its QgsGraphArcs (or edges respectively)?


texts: error while loading osm2po routing sql-script into PostgreSQL

texts:
please advice what could be wrong in my situation:

I succesfully generated sql-file using osm2po
I tried to run it within PostgreSQL and got the following error:

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

ERROR: Table 'BY2_2po_4pgr' does not occur in the search_path
SQL state: P0001
Context: SQL statement "SELECT AddGeometryColumn('','',$1,$2,$3,$4,$5, $6)"
PL/pgSQL function "addgeometrycolumn" line 5 at SQL statement

Working on Win7(x64) with PostgreSQL 9.1 and PostGIS 2.0.1
osm2po ver 4.5.25

Any ideas appreciated!


texts: How to create my own basemap in ArcGIS from JPEG files?

texts:
I have my own *.jpg tiles. How to create basemap from it and publish as service via ArcGIS for server?
I was trying to find my issue in documentation, but I didn't have success. It seems very easy for usual ArcMAP users, but not for me.


texts: Get OpenStreetMap coordinates from @id field

texts: I've downloaded some data from OSM about tourism buildings and now I need to compare them with data from another source, to check if the same elements are positioned in the same place or at least within a buffer of 10 metres around them. The other source has the coordinates of the point so it's easy, but data downloaded from osm don't have them (and also they are represented as polygons and not as points). However, in the attribute table, there's a field called "@id" which contains something like way/number. Is there a way to get point coordinates from this value in order to see corresponding points on the map and compare with the position with the other source?


texts: Problem editing raster and saving with projection and transformation intact

texts:
I am trying to make a custom tool to Arcmap using Python. It is supposed to read in a raster DEM, a point shape file and a distance, and output another raster.

The output raster should be essentially the same as the input one. In the process however, I want to edit some the raster values according to a specified function, namely those that are within the input distance from a point in the point file.

The problem with my script is the saving part. I have two different approaches. Alt 1 is working in the meaning that I get an output file. However is the projection and transformation completely off and I'm not sure how to fix it. Alt 2 is my preferred approach, but i can't figure out how to write my array to a band in the file.
Also I am not sure that the edits gets saved back to the array properly, at least they are not visible in the files I do get.

Python is still a bit new to me, so there might just be something simple that I've missed. I'm grateful for all the help i can get!

I'm using Arcmap 10.0 and Python 2.6 on Windows 7.

import os, sys, json, osgeo, arcpy, math, arcgisscripting, numpy
from osgeo import gdal, osr, ogr, gdal_array
from osgeo.gdalconst import *

inDEM   = arcpy.GetParameterAsText(0) # .tif
Outlets = arcpy.GetParameterAsText(1) # Point shapefile
outDEM  = arcpy.GetParameterAsText(2) # Path
buffr   = arcpy.GetParameterAsText(3) # Value

# Register GDAL drivers
gdal.AllRegister()

## Open the DEM
dem    = gdal.Open(inDEM, GA_ReadOnly)
band1d = dem.GetRasterBand(1)
cols   = dem.RasterXSize
rows   = dem.RasterYSize
trf    = dem.GetGeoTransform()
prj    = dem.GetProjection()
driver = dem.GetDriver()

# Get info
band_info  = arcpy.Describe(inDEM + '/Band_1')
nodata_dem = band_info.noDataValue
cellsz     = band_info.meanCellHeight
area       = int(buffr[0])
buffers    = area/cellsz # The area around each chosen cell to use for calculations

# Convert to point file to raster
outlet_raster = 'testraster'
arcpy.env.overwriteOutput = True
arcpy.env.extent = inDEM
arcpy.PointToRaster_conversion(Outlets, "Enabled", outlet_raster, "MAXIMUM", "", cellsz)
arcpy.env.overwriteOutput = False

# Convert rasters to arrays (can probably be done in a more memory-efficient way)
outl_arr = arcpy.RasterToNumPyArray(outlet_raster, "", cols, rows, None)
nodata_o = outl_arr[0,0] # I know this is not the best way, but not an issue atm
f = numpy.argwhere(outl_arr != nodata_o)

dem_arr  = arcpy.RasterToNumPyArray(inDEM, "", cols, rows, None)

# Create slope array
slope = []
x_list = numpy.arange(0,5,.01)
for i in x_list:
    slope.append(-math.exp(-i))
length = len(slope)

# Make modifications
non = []
for i in range(len(f)):
    index = f[i]
    if dem_arr[index[0], index[1]] != nodata_dem:
        r_low = index[0] - buffers
        r_up  = index[0] + buffers +1
        c_low = index[1] - buffers
        c_up  = index[1] + buffers +1
        data  = dem_arr[r_low:r_up, c_low:c_up]
        dim  = math.sqrt(numpy.size(data))
        n1 = math.floor(dim/2)
        for j1 in range(dim):
            for j2 in range(dim):
                n2 = data[j1,j2]
                dist = math.sqrt((n1-j1)**2 + (n1-j2)**2)*cellsz
                decr = length * dist / area
                if decr < length and decr != 0 and n2 != nodata_dem:
                    data[j1,j2] = n2 + slope[int(math.floor(decr))]
                elif decr == 0:
                    data[j1,j2] = n2 + slope[0]
                elif decr == length:
                   data[j1,j2] = n2 + slope[-1]
        dem_arr[r_low:r_up, c_low:c_up] = data
    else:
        non.append(f[i])

# Trying to save the changes Alt 1.
newraster = arcpy.NumPyArrayToRaster(dem_arr, "", inDEM, "", nodata_dem)
newraster.save(outDEM)

# Trying to save Alt 2.
newdem = driver.Create(outDEM, cols, rows, 1, gdal.GDT_Float32)
newdem.SetGeoTransform(trf)
newdem.SetProjection(prj)
newband = newdem.GetRasterBand(1)
newband.WriteArray(dem_arr) #


texts: Filling in An Attribute Column with Values from a List in a standalone python script

texts:
I have created an empty attribute column A via the following lines:

Layer = QgsVectorLayer("/home/usr/Desktop/File.shp", "File", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(Layer)
lProvider = Layer.dataProvider()
lProvider.addAttributes( [ QgsField("Yield",QVariant.Double) ] )

I want to fill this column with the values : t[0].Yield from my list t.
I tried by doing the following ( it does not work):

updateMap = {}
iter = lProvider.getFeatures()
idYield = lProvider.fields().indexFromName('Yield')
for feature in iter:
    Yield_fill = feature.attributes()[idYield] + t[0].Yield
    updateMap[ feature.id()] = {idYeild: Yield_fill}
lProvider.changeAttributeValues( updateMap )

Can you give me suggestions on how to edit the code above in order to accomplish what I want? I need to add information for many columns in this manner.


texts: How to convert image pixel to latitude and longitude

texts:
We are developing an android application for showing an image with geo referened with .jgw and .prj file as below. we want to identify the geo location coordinate of image coordinate(in pixel) that user touched on the image. How we can use image cordinate to find ts relative latitude and longitude using java in android.

We tried proj4 java library for android but it seems incorrect when we find the distance of two point and we couldn't trace where is the issue.

Below are our .jgw and .prj file

JGW File

0.00000180562999999
0
0
-0.00000180562999999
76.882389 (dummy location due to NDA)
8.572607 (dummy location due to NDA)

PRJ File

GEOGCS["GCS_WGS_1984",DATUM["D_WGS84",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]


texts: selecting contours based on zoom / elevation? mapserver, postgis, leaflet contours

texts:
This is what my mapserver layer in the map file looks like:

LAYER
    NAME "blah"
    CONNECTIONTYPE postgis
    CONNECTION "user=apache dbname=blah"
    PROCESSING "CLOSE_CONNECTION=DEFER"
    DATA "the_geom from (select  gid, ST_Transform(the_geom, 4326) as the_geom from public.test_blah_4326 limit 20) as subquery using unique gid using srid=4326"
    STATUS ON
    TYPE LINE
    PROJECTION
        "init=epsg:32007"
    END
    DUMP true
    METADATA
        "wms_title" "blah"
        "wms_srs" "epsg:4326 epsg:900913"
        "wms_feature_info_mime_type" "text/html"
    END
    CLASS
        NAME "platform_loc"
        STYLE
            OUTLINECOLOR 254 0 0
            WIDTH 1
        END
    END
END

The contours show up. It works. But there are so many contour lines that it just looks like a big blob on the screen.

What's the best way to "thin out" the number of contour lines based on what you're looking / how far you are zoomed in leaflet (zoom / elevation)?

I saw an attribute of the layer called PROCESSING that might help me, however, the docs aren't helpful in understanding exactly how to use it, or if it's even applicable to what I'm trying to do.


texts: How to use GRASS v.buffer function in QGIS?

texts: Trying to use the v.buffer function to create a 500-foot buffer (map units are in metres, which is why I put 152.4 meters as the buffer distance) around the boundary of a property parcel (not the centroid). I'm unsure what to populate the "'Category values' and 'WHERE conditions of SQL statement without 'where' keyboard'. All I want to do is create a 500-foot buffer around the parcel boundary, so I am finding a lot of these fields don't apply to me. When I try to run what I have entered, I get the error that there are not enough arguments.


texts:
I uploaded a zipped shapefile to Mapbox with about 500 location data points and tried to display each of those data points with one of 60 different SVG icons. Those SVG icons have been loading inconsistently on both the editor and published view.

Could this be a problem with the SVGs, the data or both?

I assigned the icons to the data points by adding a layer to the map and using the same tileset with a different filter each time. The tileset wasn't over Mapbox's capacity (it's 66KB) but maybe using it in each separate layer slows all the map loading?


texts: Network Analyst does not analyses my own created roundabouts in geofabrik shp file

texts:
I am trying to digitized the shape file data (attached) that i downloaded from geofabrik, to new roads in my region.

The network anaylst does not work beyond rounds abouts, that I create new ones. It solves the route till round about see pic 1, but beoyond that it fails. and gives me error.

Warning: No route from location "Graphic Pick 1" to location "Graphic Pick 2".
Error: No solution found.

then i removed the round about from same position, it work, so something wrong with round about, can anybody help.


texts: Oracle WMS GetFeatureInfo issue with points and lines

texts:
I am trying to do GetFeatureInfo request on MVDEMO sample data.
When it comes to polygons, it is working fine:

When I do the same type of requests on lines or points, I get an empty response:

As you can see, no attributes are returned. Any ideas why is that a case?

The attributes are enabled in MapBuilder:


texts: Time series chart CHIRPS data in Earth Engine

texts:
I am trying to derive a time series chart of the CHIRPS daily data set (precipitation) at a specific location. However, I get the following error "Error generating chart: No properties to chart."
Alternatively, it would work for me also to export the daily values in a table for the required time frame.
Please see below the script.
var lng = 69.78086;
var lat = 34.65411;
var point = ee.Geometry.Point(lat, lng).buffer(30);

var country = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
  .filter(ee.Filter.eq('country_co', 'AF'));
var aoi = country;

Map.setCenter(lng,lat, 5); // Center the map on this location, zoom level 10

var start = '2018-02-15'; // initial date of the image collection
var end = '2018-07-15'; //final date of the image collection

var dataset = ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY')
    .filterDate(start, end)
    .filterBounds(aoi);

print(ui.Chart.image.series(dataset, point, ee.Reducer.first(), 1));


texts:
Can't share my script as geoprocesing service.
The warning message: 00068 Script Script contains broken project data source: CURRENT
Here is this part of the script, where the "CURRENT" mapping module is:

# get the map document
mxd = arcpy.mapping.MapDocument("CURRENT")
# get the data frame
df = arcpy.mapping.ListDataFrames(mxd,"*")[0]
# create a new layer
newlayer = arcpy.mapping.Layer("c:/temperature/result/atlagint.shp")
# add the layer to the map at the top of the TOC in data frame 0
arcpy.mapping.AddLayer(df, newlayer,"TOP")
sourceLayer = "c:/temperature/result/symbology.lyr"
layerSymb = arcpy.mapping.Layer(sourceLayer)
updateLayer = arcpy.mapping.ListLayers(mxd, "atlagint", df)[0]
arcpy.mapping.UpdateLayer(df, updateLayer, layerSymb, "TRUE")
arcpy.RefreshTOC()


texts: What is default unit of calculated area in turf.js?

texts:
Here is a good example how to calculate area in turf.js
http://turfjs.org/static/docs/module-turf_area.html

But what is its unit of area measurement (acre/hectare) ?


texts: Cannot Edit Vertices of Polyline Feature

texts:
I am using ArcMap 10.3.1. What I have is a single polyline feature for which I am trying to edit the vertices. The only thing check in my TOC is the shapefile I'm trying to edit, so there can't be anything else selected simultaneously. I can select the feature, but when I attempt to "Edit vertices," I get an error message telling me to "Select a single editable feature to modify." I have tried stopping and restarting editing sessions, both by selecting "Start editing" under the individual entry in my TOC and from the editor toolbar, with no difference. I could just be reasonable and create an entirely new feature; I see no reason why this should not be working, and it's driving me crazy. Can anyone shed any light on what might be happening here?

The feature selected with the drop-down menu enabled.

The error I get when I select "Edit vertices."


texts:
I have a set of coordinates (points & no overlaps) in a single table,each containing attribute number fields,like in the following image.

I want to merge these points into smallest polygons by closest distance like red lines drawn around them,and aggregate attributes fields with them as well.

How can i do it in PostGIS or Spatialite ?


texts: NAD27 to NAD83, saving AFTER transformation?

texts:
I need to measure the distance (with the ruler) of points from two different layers. One layer is NAD83 and the other is NAD27. When importing the NAD27 csv file, I specified it was NAD27. Then I used: Toolbox-> Data Management -> Projections and Transformations -> Project and I projected from the input coordinate system of NAD27 to the output coordinate system of NAD83.

The issue is that I have a ton of CSV files, each with data from different projections (there is a column that specifies projection). I use Microsoft Access to connect a lot of tables and then extract the data I need, save in a table, and import to ArcGIS. The problem is that after I just transformed my NAD27 to NAD83, the lat/long in the actual table is still the same, so if I join the tables and extract data later on, I will have to do the transformations all over again.

Any advice?


texts: Stand-alone PyQGIS (QGIS 3.13) scripts with OSGeo4W and the Orfeo toolbox (OTB 7.1.0) seem to be mutually exclusive

texts:
I am changing my python environment from working with ArcGis and Arcpy to working with QGIS and PyQgis. And I also want to use the machinelearning possibilities from the orfeo toolbox. In doing so I have succesfully installed QGIS 3.13 on my Win10 PC and set up Pycharm for stand-alone PyQGIS development like explained here.
After that I succesfully installed the orfeo toolbox (7.1.0) and inside QGIS it works like a charm.
Following the instructions of the orfeo toolbox, I added the otbenv.bat to the .bat I created and the result is:
@echo off
SET OSGEO4W_ROOT=C:\OSGeo4W64

call "%~dp0\o4w_env.bat"
call "%OSGEO4W_ROOT%"\apps\grass\grass78\etc\env.bat
call qt5_env.bat
call py3_env.bat
call "%OSGEO4W_ROOT%"\OTB-7.1.0-Win64\otbenv.bat

@echo off
path %OSGEO4W_ROOT%\apps\qgis\bin;%OSGEO4W_ROOT%\apps\grass\grass78\lib;%OSGEO4W_ROOT%\apps\grass\grass78\bin;%PATH%
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT:\=/%/apps/qgis
set GDAL_FILENAME_IS_UTF8=YES
rem Set VSI cache to be used as buffer, see #6448
set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
set QT_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\qgis\qtplugins;%OSGEO4W_ROOT%\apps\qt5\plugins
SET PYCHARM="C:\Program Files\JetBrains\PyCharm Community Edition 2020.1.1\bin\pycharm64.exe"
start "PyCharm aware of Quantum GIS" /B %PYCHARM% %*

But if I try to import OTB in a stand alone scripts or aan IDE like pycharm like this:
import otbApplication as otb

the error message I receive is:
ModuleNotFoundError: No module named 'otbApplication'

The only explanation that I can think of now, is that the orfeo toolbox bindings for stand alone scripts are created for python 3.5. as is written here.
When I setup a python 3.5 environment that proves to be true. And I was able to import the module. However I could not access qgis.core from this python 3.5 environment.
And since there is no (Windows) QGIS version with python 3.5, there seems no straightforward way to combine qgis and the orfeo toolbox in stand-alone development.
I hope someone has dealt with this issue before.


texts: Wrong Data Type for my Python Script Tool?

texts:
I've got the following script and want to run it in the ModelBuilder.

import arcpy, os
dem = r'C:\folder\dem.tif'
out_folder = r'C:\outfolder'
epsgs = range(32601,32661) + range(32701,32761)
for epsg in epsgs:
    print 'Reprojecting to epsg: ',epsg
    arcpy.ProjectRaster_management(in_raster=dem, out_raster=os.path.join(out_folder,'Reprojected_epsg_{0}.tif'.format(epsg)),
                                  out_coor_system=arcpy.SpatialReference(epsg))

So I've changed the input and output this way,

import arcpy
import sys
import string
import os

dem = arcpy.GetParameterAsText(0)
out_folder = arcpy.GetParameterAsText(1)
epsgs = range(32601,32661) + range(32701,32761)

and loaded it into my toolbox

Here is the 1st problem. Did I have chosen the wrong Data Type?
I want the script to save the created files into the existing folder.

Then I tried to write another folder to the path and the error disappeared.

Y:\Home\hoffmada\World_rivers\Verarbeitete_Daten\Flowdirection_UTM_120\flowdirection

But when I want to execute the script I got this error message.

So what Data Type do I have to choose to write down the new rasters into my folder I've chosen?


texts: Create new layer and populate attribute table from QtWidget PyQGIS

texts:
in my QGIS plugin, one of the dialog has a QtWidget table where the user can add items to the columns.

Each column of the table is a single list and the final result is a dictionary where the keys are the column name and the values are the single entries, like:

# example of single lists
sp = [1, 2, 3]
th = [0.2, 0.3, 0.1]

#example of dictionary
dic = {'stress': sp, 'thick': th}

From the dictionary I need to create a memory layer (no needs to have a spatial layer, just a table one).
This memory "table" should have an attribute table where the files are the dictionary keys and the features are the dictionary values:

====== =====
stress thick
====== =====
1       0.2
2       0.3
3       0.1

Any suggestion?

EDIT

I made some progress.

from PyQt4.QtCore import QVariant
# create a memory layer
vl = QgsVectorLayer('Point', 'layer' , "memory")
# data provider
pr = vl.dataProvider()

# lists of values, dictionary and list of field type, just as example
sp = [1, 2, 3]
tt = [10, 20, 30]

dd = {'stress':sp, 'thick':tt}
var = [QVariant.Int, QVariant.Int]

# create the fields of the layer with the dictionary keys
for i in dd.keys():
    for j in var:
        pr.addAttributes([QgsField(i, j)])

# update the fields
vl.updateFields()

until here is fine, but what I'm not able to do is to loop over the dictionary values and add them to the attribute table.

The workflow should be like this:

# start editing
vl.startEditing(

ft = QgsFeature()

# new attributes, if just one for each field
ft.setAttributes([1,10])

# send to data provider
pr.addFeatures([ft])

# commit the changes
vl.commitChanges()

But, how can I add many features in a single setAttributes method?


texts: Using render_list to render only map extract

texts:
I am trying to render an area where I have added an extract from the planet file with render_list. The command I am using is

render_list -x 159 -X 179 -y 160 -Y 410 -n 2 --socket=/var/run/renderd/renderd.sock --min-zoom=10 --max-zoom=10

to render zoom level 10. I already rendered 0-9 zoom levels for the whole planet, but want to render only the areas where I imported maps for the bigger zoom levels.

Running the command above gives message "Starting 2 rendering threads" and stays for a long time without doing anything. Any ideas why?


texts:
I´m looking for blogs or good websites about digital/web cartography...i´m not really interested in the deep technical ramblings (i´ve the forum lists for that) but rather inspirational blogs which showcase great digital cartographic techniques or projects and how they are created.

If there are such things here in the web, please let me know...I´,m also interested in statistical visualisation blogs like these..http://datavisualization.ch/showcases/, http://flowingdata.com/...but I´ve never seen such blogs which specialise in digital/web cartography.


texts: Geoserver layer from the ArcGIS image mosaic

texts:
I've received a number of GeoTiff files most likely created in the ArcGIS in the EPSG:32641. Here are *.tif files itself, *.aux.xml, *.tfw, *.ovr.

I want to put them to the GeoServer as image mosaic data source and serve the applications with the layer through WMS. So I've just added the image mosaic data source and pointed the files directory.

Here is an issue with creating of the ImageMosaic layer in the GeoServer: when trying to preview it in OpenLayers I can only see picture/colors only for several zoom levels and only some tiles and only after clicking. But most of the image is white. However if click to the any white point of the image its possible to observe valid RGB band values below the image:

The actual data size covers the city. However data looks fine in the ArcGIS Explorer Desktop.

Could it be possible that provided data is invalid or wrongly formed and inappropriate for GeoServer? Or am I doing something wrong?

I've also tried to join the .tif and all other related stuff (.tfw and so on) and add the result tiffs as single GeoTiff images as well as warp to the EPSG:4326 but the result is same.

Probably it is same issue as described at link but I'm not sure.

UPD: I've tried to add my geoserver WMS layers as a source GIS server into ArcGIS Explorer Desktop. The result is following:

So looks like on the different zoom levels the color data is visible for some areas. For other areas on the same zoom only contours are visible.

UPD2: Merged to one big (120mb) GeoTiff the output looks better but on the closest zoom the problem still exsits.


texts: How do I get accurate DTM from DSM to get tree height?

texts:
I have drone imagery which was generated in Agisoft. From the imagery a DEM was generated and exported.

I would like to know how do I get the Digital Terrain Model from this DEM? as I currently have the Digital Surface Model?

The method I have tried doing was creating points at different locations in between the trees(on the ground) and interpolating the points to create a DTM by subtracting the DSM and DTM (DSM-DTM= tree height) from each other.

However the results are not consistent enough to achieve tree heights accurately.

I do have ArcGIS Pro.


texts:
I want to generate WFS and WMS from a PostGIS database. To do this I want to make a script with Python (or with another language) to store my shapefile on GeoServer and then deploy services automatically.

Is there any good approach to follow?

OWSlib could be one?

And what about geoscripting?

Point of situation:
I have made a script for creating shapefiles and generate metadata so the next step is to build a script that post my datasets into GeoServer to generate services. The number of datasets are constantly growing because is periodic data, and I will assume that this upload script simply will be executed occasionally by user interaction or eventually a daemon could be configured.

The data will be provided from different database, converted into specific datasets(shapefile) and will feed a GeoServer (or MapServer), service.

Is there any script that already make this interaction, to upload shapefiles from a specific directory into GeoServer?


texts: Detecting changes between two shapefiles on basis of attributes as well as spatial location in single query?

texts:
In MapInfo we can easily compare two files on the basis of sql query where map info works on rows. You need to update only latitude and longitude and create points on the basis of lat/longs in both the two files.

Now run a sql query in MapInfo -

Where Fileone.obj intersect file2.obj and file1.floodattribute1 = file2.floodattribute2 = Output

will give you all the matching records just invert selection and you will get only change areas either on spatial or attribute level.

Can similar activity done in ArcGIS?


texts: QGis 3.4: GPKG vector layer with object form (relations with spatial and non-spatial datasets) to postgres database

texts:
I have QGis project with a geopackage included. One layer in the geopackage contains a attribute form with relations to other layers (spatial and non-spatial). I have to convert it to a postgres database. Is there a way to transform the attribute form (as much as I know created with drag and drop) to a postgres db? And then to import the data?

I can create the database by transforming the geopackage to a postgres with GDAL ogr2ogr but how can I automatically include the relations in the attribute form?


texts:
So, I don't really have much knowledge about the coordinate system and sometimes I get stuck in certain things that I just can't get around it. For example, whenever I'm working with a vector that is in UTM and I need to see a raster of that place, I use the plugin to see the google satellite in the background.
But how can it be perfectly placed with the vector? Wasn't it supposed to have like a difference between the two since the vector is in UTM and Google works in another coordinate system? How does that work? I am doing something wrong?
Could someone explain that to me or give me some references that I could be reading to fully understand the subject?


texts: Why can't lines plotted on google maps suffice for legal determination of property boundaries?

texts:
If there is a way to calculate a maximum error in a line plotted on google maps (road map, or satellite view) , then that should be admissible in a court of law, shouldn't it ?

For example, my deed describes my property boundaries starting from a "point of commencement" and then gives distances and angular measures to eventually arrive at the POB and similarly then continues to outline the actual boundaries.

If I locate the "point of commencement" I could, with a bit of coding, trace exactly those lines on google maps.

Is there some way to estimate the maximum error in real world distances of such lines as shown on the map from their actual location?


texts: Check direction of linestring against corresponding polygon using SpatiaLite

texts:
I have a polygon and a linestring that is part of the polygons boundary, but not necessarily digitized in the same direction as shown here:

In this case the linestring runs in the opposite direction of the polygons boundary.
How can I check if the linestring has the same direction or not using SpatiaLite?


texts: Plot RGB values for each feature seperately

texts:
I want to plot each rgb value for its respective feature of the layer via pyqgis. Atm I struggle to properly implement it for each layer individually. With my current code below, all features get the last color code in row no.3. and therefore they all have the same color instead of every element ist individual one.
I use QGIS 2.8.8 Wien and therefore proposed solutions for QGIS 3 do not work properly.
b = iface.activeLayer()

symbols = b.rendererV2().symbols()
symbol = symbols[0]

for f in b.getFeatures():
    expr = f['rgb']
    express = "color_rgb("+str(expr) +")"
    symbol.symbolLayer(0).setDataDefinedProperty('color', express)

b.triggerRepaint()
iface.legendInterface().refreshLayerSymbology(b)


texts:
This is a follow up to Output field name not showing up in spatial join result in which I received no response.

I am trying to convert a shapefile to a feature class, I also have a GeoJSON file with the field names as I would like them to appear in the output feature class. The arcpy.FeatureClassToFeatureClass_conversion functionality accepts a field map, so I tried to build one by looping through a GeoPandas geodataframe of the GeoJSON file (to set my output field names) and by adding the shapefile to a fieldMappings object, and replacing the field map for each iteration:

import archook
archook.get_arcpy()
import arcpy
import geopandas

file_path = 'damage_results/Culverts/Culverts.shp'
original_fields = arcpy.ListFields(file_path)
original_fields = original_fields[2:]
gdf = geopandas.read_file('damage_results/Culverts/result_culverts.json')
new_field_names = gdf.columns[:-1].to_list()

fms = arcpy.FieldMappings()
fms.addTable(file_path)

for i,field in enumerate(original_fields):
    print i
    fm = fms.getFieldMap(i)

    out_field = fm.outputField
    out_field.name = new_field_names[i]
    fm.outputField = out_field
    fms.replaceFieldMap(i, fm)

for fm in fms:
    print fm.outputField.name
arcpy.FeatureClassToFeatureClass_conversion(file_path, 'intermediate_damage.gdb',
                                            'Culverts', field_mapping = fms)

The second for loop is to verify that the output field names have been written properly, and it seems from within the script that they are. However when I open the Culverts feature class the field names are unchanged from the shapefile.

Does the output field attribute in the FieldMap object work?


texts:
I try to seed a wms service that is very very slow. Therefore the seeding fails with the message "Operation timed out after 600000 milliseconds with 0 bytes received".

It seems not to be possible to set the time out time. I found the connection_timeout option in the wms-source of the mapcache.xml-file but it not work and changed the value from 600 to 1800:

<connection_timeout>600</connection_timeout>

to

<connection_timeout>1800</connection_timeout>

but the error message remains the same "... 600000 milliseconds ...". The change seems not to be adopted. Even after a restart of the apache webserver.

The command for the seeding

mapcache_seed  -c mapcache.xml -t osm_eu_3035 -e 542100,-369400,7530900,6619400 -z 0,10

The osm_eu_3035 service uses the local very slow wms service.

What did I wrong? Is there any other possibility to increase the time out time?


texts: Python - OGR: Distance different than in ArcGIS

texts:
I want to determine the shortest distance between a point and a line. I use the following Python script to do so. If I compare the result of the script with the result of ArcMap's Near Analysis, I get a completely different results. The shapefiles are in the CRS GCS_North_American_1983 and projected with Albers. What could be the reason for the difference?

import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')

lineshp = driver.Open('U:/My Documents/Tool/shp/line.shp', 0)
linelyr = lineshp.GetLayer()

pointshp = driver.Open('U:/My Documents/Tool/shp/point.shp', 0)
pointlyr = pointshp.GetLayer()

line = linelyr.GetNextFeature()
point = pointlyr.GetNextFeature()

line_geom = line.GetGeometryRef()
point_geom = point.GetGeometryRef()

dist = point_geom.Distance(line_geom)

print dist


texts: I'm trying to display some weather data in ArcMap 10.4.1. Starting with a blank map, I've successfully displayed the data from a GRIB file. I've figured out the lat/lon grids and those are showing. How do I add country borders to the map?


texts: ArcGIS Pro pixellates vector marker symbol when published as web feature layer

texts:
I've added multiple SVG images to a feature layer as custom symbols in ArcGIS Pro. While these show up exactly as expected in ArcGIS Pro, the same symbols appear pixellated once the feature layer is published to ArcGIS Online.

I am inexperienced in both designing SVG images and creating custom symbols in ArcGIS so I do not understand whether there is a problem with the ArcGIS settings or the way I have created the SVG file.
Is it possible to publish custom vector symbols within feature layers to ArcGIS Online?

I managed to get some slightly better results by adding the custom vector symbol as a new style in ArcGIS Pro - that is, ArcGIS Pro seems to display the symbol in its own menus in a way that looks more like one of the default symbols. However, whenever the feature layer is published online the symbol is still pixellated. Seems like a low-quality conversion from vector to raster during the upload process.


texts: Isolating specific line of multiline feature class in ArcGIS Pro

texts: I have a polyline feature class containing conjoined line features. Each line is composed of multiple individual lines separated by vertices. I would like to either isolate the central line in the image and extract it, or remove the tributary lines that are running perpendicular to the centre line as in the image. I have no attributes to select on but do have a bounding polygon that could potentially be used. The yellow dots in the image represent the start and end vertices of each line. Manually, I can select each individual line that I don't want and delete but the eventual aim of this is to automate within a Python script. The 'Trim line' tool in ArcGIS Pro does what I want but I have to run it tens of times and also trims the ends of the central line that I want to keep.


texts: I am working in QGIS 4.8, and I am trying to use the Whitebox Processing plugin. I have Whitebox Tools and Whitebox GAT zip files installed and unzipped on my drive. However, I cannot use the Whitebox Processing plugin, and Whitebox isn't showing up as a provider. Has anyone else faced this issue?


texts: Sharing QGIS data with ArcGIS Desktop user?

texts:
I mainly use QGIS and have a large project to complete. Once I am done I'll need submit my data to an ArcGIS user. The data will consist of mostly point and polygon shape files.

How do I set up a QGIS project from the start if I know that I'll be handing it over to someone who uses ArcGIS Desktop?


texts:
I've just added a couple of more WFS backed vector layers to a map and the performance has dropped.  Whilst debugging this I've noticed that openlayers is considering every feature in the visible layers when rendering.

If I'm zoomed right in so there are a only couple of features from a few layers and then say pan the map up a small amount.  I'd expect openlayers to use a spatial index to identify which features need to be rendered and only render those.

However, I'm seeing my context based styles being asked for symbolizer info for every feature in every layer.  E.g. my externalGraphic function is been called for features that are on the other side of the country!

  externalGraphic: function(f) {
      return map.zoom < switchZoom ? icon : "/grails-app/feature/image/" + f.fid;
  }

When zoomed out my layers are not visible but when zoomed in there are suddenly 1000s of feature being considered.  Is there a way to get openlayers to use some filter when rendering?  I had thought about adding map zoom/pan event handlers to map features in and out of layers using a spatial filter but I'm pretty sure I shouldn't have to do that.

I don't believe that the fact that my data originally comes from a WFS source is the issue here.  My app also runs offline, the data, once retrieved is stored locally, the retrieval and resources needed to hold these in memory is fine.

There must be a way of getting OL to apply some intelligence when rendering vector features when zoomed in, it shouldn't be even considering rendering features that are hundreds of miles away from the current bounds, this is triggering retrieval of images from either the local file system or the server when I the zoom level switches from a generic icon to a more detailed device specific image.

In total I have about 5000 features, this should be a trivial amount.  I've worked around this by handling move/zoom events and mapping my features in and out of the layers using a quad tree from the bespoke library I wrote years ago and want to move away from.  This is an ugly hack that I find it hard to believe that that there is not a more elegant way to handle, I'd expect to provide the map with my vectors and let it manage them.


texts: Inquiry of using sDNA in analysis of park accessibility in an island and its accuracy

texts:
I have a question  whether sDNA closeness analysis could be used in any region/ city. I am working on an analysis of  the accessibility of parks in Manhattan. You know, Manhattan is a narrow linear island so the core area with high integrity (closeness) sometimes is pretty close to the boundary of the island with low integrity.  It resulted in low accessibility of the parks in the edge of Manhattan when performing calculation with the sDNA model. However, the pedestrian flows are high in parks along boundaries according to observation data. So I am wondering if the methodology is more appropriate to use in continental rather than island sites.
Have you done any research about this limitation?
Attached are the diagrams I drew with sDNA and GIS. They are local closeness at R 800m, accessibility of park at R 800m (join with the closeness of road). The other diagram with red dots shows the times of crowd observed at those parks. You can tell many dots distributed in the parks along the boundaries, which can't be explained by the closeness of the road system.
I am wondering if closeness/integrity is a suitable  method to explain the crowd data.
Update: The closeness I mean in the analysis is actually NQPDA, which is the model you said better to describe accessibility. So the first and second suggestions are not working. Adding scenic weighting/ program weighting may work since their distribution reflects such apparent patterns.
How to add weighting in integral analysis in sDNA?
I think I may use the "learn" and "predict" function in sDNA.  I plan to use three weighting factors (NQPDA, population density nearby and number of programs) to run a predicted model for estimating the possibility of crowds.
How to apply these three  factors in "learn" and "predict" function in sDNA?
I looked at the sDNA manual but it didn't introduce details. https://sdna.cardiff.ac.uk/sdna/wp-content/downloads/documentation/manual/sDNA_manual_v3_0_alpha1/guide_to_individual_tools.html#learn


texts: Cost Distance Analysis QGIS with Vector files

texts: I am trying to calculate the costpath from a location (polygon) to a point dataset using QGIS. I found two tools (r.cost.raster & accumulated cost (isotropic)), which both would work, if I would use only raster datasets. Unfortunatly one of my Datasets is a polygon so I can't use these tools. Any ideas how I could work around this problem?


texts: Open Street Map and PostgreSQL what do I exactly need?

texts:
I make my first steps in the GIS world and everything is very new for me even if I work on a "huge" project (ie: for a single programmer) very related to maps.

My project needs to know informations below for any location of the world:

Road type
Road name
Road coordinates
City boundaries
Water way name
Water way coordinates

I found that I could have all this informations thanks to Open Street Map and my project begins to work thanks to this. But now I need to work on bigger maps (ie: bigger than the 50.000 XML nodes limit OSM export provides). I know I need to put the whole planet in a database.

There is a lot of documentation about PostgreSQL and Open Street Map on the Internet but I don't know where to begin.
In fact I found documentation about Osmosis in order to put a planet dump in a database, but I also found about osm2pgsql. The point is I don't know what is the best for me.

It seems (but I can be wrong) osm2pgsql has been made in order to render mapnik maps, thing I don't need at all because I have to know roads/cities/water ways informations... But I am not sure to have correctly understood what mapnik can do.

My question is: What can I use to insert OSM datas in a database in order to reuse it later as easily as an OSM exported file?

Thank you.

EDIT: The computer I use is a dedicated server (Intel(R) Pentium(R) E2180 2.00GHz and 2GB of RAM) running a GNU/Linux Debian.

PS: I know my question seems dumb, because osmosis can insert but also extract a xml file from a database. Anyway, I need to know how you would do it.

PS: If someone knows where to find nice OSM+database tutorials for beginner I would appreciate!


texts: Does anyone know if it's possible to generate line (the distance between each one is 2 meters) in a polygon?


texts:
I want to compare road centreline attributes to see if they have different attributes where they meet.  For example, perhaps the split in road lines is because the speed limit changes. In this case I want to compare the two lines to test whether the speed limit changes at this split.

My initial idea was to use LineJoiner to merge connected lines, using the speed limit attribute to do the join.  Then to create points at the start and end points for each line, and then to compare the two closest points to see if the speed limit changes.  However, I am struggling in creating the points for the start/end locations and keeping the attributes.

The vertexExtractor only seems to create one point at the split locations, which does not allow the comparison to be done.

Are there better transformers to use to achieve what I am trying to do?


texts: How to remove zonal anisotropy from directional variogram using gstat in R?

texts:
I have a dataframe named seoul315 contains the 107 data of PM10 concentration at 1 march,2012,5.00 am. please, download. I have tried to plot semi-variogram for kriging bt the following code:

library(sp)
library(gstat)
library(rgdal)

seoul315<-read.csv("seoul315.csv")
seoul315<-na.omit(seoul315)

hist(seoul315$PM10)   #Data are not in Gaussian form.So need to log transform
hist(log(seoul315$PM10))

coordinates(seoul315)<-~LON+LAT
proj4string(seoul315) =  "+proj=longlat +datum=WGS84"
seoul315<-spTransform(seoul315, CRS("+proj=utm +north +zone=52 +datum=WGS84"))

#plot Omnidirectional Variogram
seoul315.var<-variogram(log(PM10)~1,data=seoul315,cutoff=80000, width=10000)
seoul315.var
plot(seoul315.var, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 315")

#Model fit
model.315<- fit.variogram(seoul315.var,vgm(psill=0.11,model="Gau",range=60000/sqrt(3),nugget=0.04),
                          fit.method = 6)
plot(seoul315.var,model=model.315, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 315")

#Directional Variogram
seoul315.var1<-variogram(log(PM10)~1,data=seoul315,width=10000,cutoff=80000,
                         alpha=seq(0,135,45),tol.hor=15)
plot(seoul315.var1,model=model.315, cex=1.1,pch=16,col=1,
     main="Anisotropic Variogram for PM10")

I got the directional variogram like this:

But from this variogram I can see sill are not same in all direction that means its zonal anisotropy. I think in this case I can't use anis() function which I used to remove the geometric anisotropy. How can I remove the zonal anisotropy for using the variogram in kriging? Is it possible by using gstat in R?

Actually I can't find the proper way for doing this after rigorous searching. If you need further information please let me know.


texts: Make Level-2 Landsat images comparable with Level-2A Sentinel-2 images

texts:
Is there a way to make Level-2 Landsat images (e.g. Landsat-5 images ordered from EarthExplorer) comparable with Level-2A Sentinel-2 images?

The objective is to calculate water or vegetation indices (NDVI or NDWI) from both and be able to accurately evaluate the change between the two, since they were not acquired on the same year obviously. The Sentinel-2 Level-2A images were either downloaded from Sentinel Hub or atmospherically corrected with Sen2cor.

While both of them are atmospherically corrected, the pixels values I find in their image bands are on a completely different order of magnitude.


texts: Error with LinearFit in geographical line

texts:
I have a set of pixels that follow a straight line.  However, the values aren't all clearly visible.  I want to create a Line Feature that I can later use to identify it.

I have created a Lon Lat image and have masked it with the value of my pixels so that only the Pixels of interest have a value.  It should be quite simple, having 'longitude' and 'latitude' as bands to create a linear fit of these geographical points.  firstLine

However, when I try to do a linear fit of this cloud of points I receive an error.  (Code is below)

Image.reduce: ReduceBands reducer must have one input: LINEAR_FIT.

var lonLat = ee.Image.pixelLonLat();

var myMask = myPixelLine;

var firstLine = lonLat.updateMask(myMask);

var linearFit = firstLine.select(['longitude','latitude'])
                            .reduce(ee.Reducer.linearFit());


texts:
I am attempting to package QGIS from source and the creatensis.pl script is exiting with the error:

  Running NSIS
  Can't open script "Master'"
  running nsis failed [makensis -V1 -DVERSION_NAME='Master' -DVERSION_NUMBER='2.9.
  0' -DBINARY_REVISION=1 -DVERSION_INT='2090001' -DQGIS_BASE='QGIS Master' -DINSTA
  LLER_NAME='QGIS-OSGeo4W-2.9.0-1-Setup-x86_64.exe' -DDISPLAYED_NAME='QGIS 'Master
  ' (2.9.0)' -DSHORTNAME='qgis' -DINSTALLER_TYPE=OSGeo4W -DPACKAGE_FOLDER=osgeo4w/
  unpacked-x86_64 -DLICENSE_FILE='osgeo4w/unpacked-x86_64/apps/qgis/doc/LICENSE' -
  DARCH='x86_64' QGIS-Installer.nsi] at creatensis.pl line 431

UPDATE: Looks like the creatensis.pl is expecting a NSIS script called Master.

Why is makensis expecting looking for a script called Master?  The command looks like it it should be looking for the QGIS-Installer.nsi script.

I am using Windows 8.1 64 bit.


texts:
I have a number of tables that contain known values in a field. For example, I have a soil type table with one field in it. That field contains a list of known soil types. I have another table that has a list of known vegetation types. I have various other tables that contain lists of other known types of other data.  What I want to do is:

Read in each table, do a cursor search and add the values to a list
via tableList.append
Now that I have built each list, I want to compare the values from a
separate feature class table to the known values in the list I have just created using if row.SOIL_TYPE in SoiltableList:
In order to do this, I need to reset the workspace I guess, because
you cannot set a path parameter for listing data in any of the arcpy
Listing methods.

env.workspace = r"C:\BASE_TABLES.gdb"

SoiltableList = []
VegtableList = []

for table in arcpy.ListTables():
   if table == "SOIL_TYPE":
      for row in arcpy.SearchCursor(table)
          SoiltableList.append(row.SOIL_TYPE)
   if table == "VEG_TYPE":
      for row in arcpy.SearchCursor(table)
          VegtableList.append(row.VEG_TYPE)

# more iterations of table

# Reset workspace
env.workspace = r"C:\TEST_GDB.gdb"

for fc in arcpy.ListFeatureClasses():
   if fc == "SOIL_FEATURECLASS":
      for row in arcpy.SearchCursor(fc):
         if row.SOIL_TYPE in SoiltableList:
            print "YES"
         else:
            print "No"
  # More iterations of if/else logic to compare values from fields from various feature classes to the known values in the list I built.

The 2 questions I have are:

Do I need to reset the workspace the way I did above? is there a better way of doing this or is what I am doing completely acceptable?
Does my code make sense? Am I over complicating what I am trying to do? I though that maybe I would define a function or 2 to simplify the code, but I'm still in the process of figuring out what else I want the code to do.

Thanks for any suggestions.

Mike


texts: I am working on a project where I need to find the mileage between an inspector's daily stops. I have a route of 11 stops and am having trouble figuring out if its possible to export out the mileage between stops but not the cumulative mileage


texts: How to fill the color in choropleth map in absolute value?

texts:
Problem

I am planning to publish an article about 'how many votes are counted in each region' every one hour after the June election in Korea is over at 6 p.m. So at 7 p.m., at 8 p.m., at 9 p.m., and so on, I am going to publish an article with this infogram below. Plus, the stronger the color, the bigger percentage of total poll is counted.

So in my code below, the column 'ESRI_PK' will change every one hour, to be more precise will go up, since people will count more votes. (here, ESRI_PK plays the role of the 'vote counting rate', even though in reality it isn't.)

What I would like to do is to show how many votes are counted cumulatively as time passes by. However, the scenario on the left hand side('As-is') illustrates what my current code is doing. Even though the cumulative vote counting rate goes up, the region A has always the lightest color because out of the 3 regions, the A region has the lowest counting rate.

I wish my code would follow the scenario on the right hand side('To-be')! Unfortunately though, the current plotting method will only do the left hand side.

My code

import geopandas as gpd
import matplotlib.pyplot as plt

plt.xkcd()
seoul=gpd.read_file('../json/seoul_municipalities_geo.json')
# You can download the same file from the Github below
# https://github.com/southkorea/seoul-maps/tree/master/juso/2015/json

final_pic=data_result.plot(figsize=(14,10),linewidth=0.25, edgecolor='black', column='ESRI_PK',cmap='Blues',scheme='quantiles',legend=True)
for index,row in seoul.iterrows():
    xy=row['geometry'].centroid.coords[:]
    xytext=row['geometry'].centroid.coords[:]
    plt.annotate(row['SIG_ENG_NM'],xy=xy[0], xytext=xytext[0],  horizontalalignment='center',verticalalignment='center')
    plt.axis('off')


texts: Clustering points/polygons based on proximity (within specifed distance) using QGIS?

texts:
I have seen similar questions asked about grouping points and Concave Hull being given as an answer, but in my case there are no "features"/attributes that will classify the points. The only parameter I want to group the points is through proximity to each other.

I have over 800,000 points which are spread across with some close to each other. I want to create polygons that cover points that are within 500m from each other, and lone ones will have a polygon of their own.

Clarification: The clusters need to have a maximum specified radius. Clusters themselves should not merge based on how far apart they are to each other.

Could you please let me know how I can approach this?


texts: standalone python script unable to execute merge vector layers

texts:
I am trying to run the qgis merge vector layers tool in a standalone script however I am getting the following error:

Unable to execute algorithm
Wrong parameter value:/Users/jt/workspace/Manna/workflow/mesh.shp;/Users/jt/workspace/Manna/utm/CW_roads_clean.shp

I first ran the tool normally and then copied the code out of the history. I have also successfully ran the code in the qgis python console. It is just when I try to run it in the following stand alone script that it has issues.

import sys
sys.path.append('/Applications/QGIS.app/Contents/Resources/python/plugins')

from qgis.core import *
from PyQt4.QtGui import *
import processing
from processing.core.Processing import Processing

app = QApplication([])
QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()
Processing.initialize()
Processing.updateAlgsList()

processing.runalg("qgis:mergevectorlayers",
    "/Users/jt/workspace/Manna/workflow/mesh.shp;/Users/jt/workspace/Manna/utm/CW_roads_clean.shp",
    "/Users/jt/workspace/Manna/ten_meter/temp_merge.shp")


texts: Converting ASCII data to GeoTIFF for ArcGIS Desktop?

texts:
I have a data file that gives me a porosity value and the lat/long that it applies to as a box.

For example

Porosity=14.58
90.0    0.0
90.0    360.0
89.0    360.0
89.0    0.0

I need this data in a format that ArcGIS can read and display as a map, preferably as a geotiff. What do I need to do to put this in a format that Arc can understand?

I have access to ArcGIS 10.4, Python, and MATLAB.


texts: Clipping a raster DEM between two elevations

texts:
I am trying to clip a raster DEM between two elevations, specifically a portion of the contiguous United States in an SRTM 90m DEM.

I would like to create a vector polygon mask that overlies the all parts of the DEM that lie between two elevations; for example, between 1000m and 2000m.

At the moment I have downloaded the SRTM section I need in both ASCII and GeoTIFF formats. I am using QGIS 2.18.15.


texts:
I'm trying to figure out an efficient way of pulling data about ground cover along a specific path, where I know the starting lat/long and the end lat/long.

I've found that the USGS has ground cover data available in GeoTIFF files, but I'm having a hard time figuring out how I can get data from these files along a path easily. Before I go down the long road of learning about this format and trying to write something to efficiently get this data, is there a better way? Am I making a lot more work for myself?

I already get the topography along the path by using the Bing elevation API, but there doesn't seem to be any API out there I can find for ground cover. Specifically, I want to supplement the topography with a 'best guess' about tree coverage, etc, along the 2D path that I'm drawing with the Bing data.

I'm currently using a PHP backend and displaying the data using d3.js on the client side (this is for a web application) but I don't really care how the data is found, as long as I can do it in an automated way. So I'd prefer a command line tool or some kind of library for PHP/Python/Ruby/etc.


texts: ArcGIS JavaScript runtime map doesn't display

texts:
I started following the code sample the Esri provided, but when I run it the map doesn't show up.

This is the link:
https://developers.arcgis.com/javascript/latest/sample-code/get-started-mapview/index.html

This is the code:

<!DOCTYPE html>

<meta charset="utf-8">
<meta name="viewport" content="initial-scale=1, maximum-scale=1, user-scalable=no">
<title>Get started with MapView - Create a 2D map</title>

<link rel="stylesheet" href="https://js.arcgis.com/4.1/esri/css/main.css">


texts: Are there any other providers of map data in the UK apart from the Ordnance Survey?

texts: The Ordnance Survey charges a lot for most data, so what are the other options?


texts: Calculating fields (multiple) with different conditions for each field using ArcGIS Pro

texts:
I have a population dataset where I want to reclassify the age groups and include gender. I have the field "gender", 1 for men and 2 for women, and the field "age". One row for each person in the dataset.

I want 5-year age groups based on gender: Men 0-4 years (M_0_4), men 5-9 (M_0_9)...up to M95_99 and W_0_4, W_5_9 and so on. I can do this manually for each field with this code:

M_0_4=
    Reclass(!gender!, !age!)

    #codeblock:
    def Reclass(gender, age):
       if (gender = 1 and age >= 0 and age <= 4):
            return 1
        else:
            return 0

Is there anyway to do this for all fields using calculate field (multiple) in one go? Each field would need a different condition. I have tried making a separate expression for each field and then just copying the codeblock for each and changing the expression name and conditions without luck, e.g:

#expressions:
M_0_4=
  m_0_4(!gender!, !age!)
M_5_9=
  m_5_9(!gender!, !age!)

#codeblock:
def m_0_4(gender, age):
   if (gender = 1 and age >= 0 and age <= 4):
        return 1
   else:
        return 0

def m_5_9(gender, age):
   if (gender = 1 and age >= 5 and age <= 9):
       return 1
   else:
       return 0


texts:
I've got a function that takes an Openlayers mouse click and puts the click location into a SQL statement which then gets passed to a JSP script which runs the SQL statement on the PostGIS database and returns geojson. I've grown bleary eyed googling for possible reasons why the line var features = geojson.read(e.responseText) always returns features as null. I'm interested in getting an attribute/property of one of the features in the returned geojson.  What am I not understanding about reading the returned geojson?

function getWatershed(e) {
    var geojson = new OpenLayers.Format.GeoJSON();
    // Get the mercator coordinates from the viewport
    var xy = map.getLonLatFromViewPortPx(e.xy);
        // Find the watershed ID.
        var sql  = "select st_asgeojson(the_geom), id, area_sqmi ";
            sql += "from medford.subshed2 ";
            sql += "where st_contains(the_geom,st_transform(st_setsrid(";
            sql += "st_makepoint(" + xy.lon + "," + xy.lat + ")";
            sql += ",900913),26918))";
        OpenLayers.Request.GET({
            async: false,
            url: "./04-sql-to-json.jsp",
            params: {sql: sql},
            success: function(e){
                var features = geojson.read(e.responseText);
                if(features) {
                    if(features.constructor != Array) {
                        features = [features];
                    }
                for(var i=0; i<features.length; ++i) {
                    watid = features[0].attributes.id;
                }
                }
            }
        });

The geojson returned from the server looks like this (the coordinates were abbreviated):

{"type":"FeatureCollection",
"features":[{"type":"Feature","geometry":{"type":"MultiPolygon","coordinates":[[[[251083.32142683025,4381268.9660757259,0],[251139.94862688333,4381204.021947667,0]]]]},
"properties":{"area_sqmi":"344.0191",
                     "id":"PU2_3900_3750"
} }
]}

What am I not understanding about reading the returned geojson? Why is features always null?


texts: How is the default label field chosen by ArcMap?

texts:
I have 2 questions:

I know how to manually go to Layer Properties>Labels tab>Label Field and choose my label field.  I also know that I could set up my primary label field and save the setting as a layer (*.lyr) file.  This question, though, is more about how the label field is chosen by the program prior to user interaction.

First: When a new layer (shapefile, ArcSDE layer, personal/file geodatabase featureclass) is added to ArcMap, how is the default label field chosen?  From my experience, it isn't chosen by field order, nor alphabetical by field.  My best guess so far is that it might chose the first field of text/varchar (alphanumeric) type.

Second: Is there a way to define what the "primary label field" will be for a layer within the shapefile, sde, or gdb files, or are you stuck with saving layer files or letting the system pick for you.


texts: Converting ArcGIS topology rules file (*.rul) to *.txt, *.csv or *.xls for reviewing?

texts:
Has anyone converted ESRI's topology rules file (.rul) into another format?

I built a complex topology rules set and would like to QC it in something easier to read and work with than ArcCatalog.

A text file and spreadsheet format would work great.


texts:
I cannot understand what the Value Relation widget of the Fields menu in QGIS actually does.

For example: I have a point layer with some columns (ID, value1...) and a table (so just a dbf file) with other values (ID, valueA,...). If I open the Fields menu of the point layer I can set some options for the selected column: Layer (set to table.dbf), Key column (set to ID) and Value column (set to valueA).

Now, if I create a new point, it appears in the row of the widget column a dropdown menu where I can choose the valueA of the table.dbf. And that's ok.

But, what is the meaning of the Key column?

Does someone have an explanation for that? Am I missing something?


texts: ArcGIS map server not loading layer with ArcGIS API for JavaScript?

texts:
I am trying to load the regions from here http://apps.geoportal.icimod.org/ArcGIS/rest/services/Nepal/VDC/MapServer,
but nothing is showing up.

Is there something wrong with my code?

<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="initial-scale=1,maximum-scale=1,userscalable=no">
  <title>MapImageLayer - 4.1</title>
  <link rel="stylesheet" href="https://js.arcgis.com/4.1/esri/css/main.css">


texts: Mathematical rule to convert between degrees and centimeters based upon projection

texts:
I have some Tanzania drone imagery data in a local geographic projection EPSG:32737 that uses the unit of meters. Then I have other vector data that is in EPSG:4326 which is denominated in degrees. In order to map the vector data on the imagery, I converted the imagery over to 4326. Since I believe 4326 is the most popular projection, I figured having all data in that projection would make things easier and avoid inadvertent projection errors.

The problem is that I want to resample the imagery to a lower resolution. In the original projection the units were meters but in 4326 the units are degrees. So my question is, is there a simple rule of thumb to convert between degrees and meters--given some local projection? Like what is the algorithm to compute this. Note the total study area is about 2 miles x 2 miles in size, since the resolution is so high.

Note, I understand that a degree corresponds to a different number of meters, depending on where on earth the imagery or geography is captured. So that is fine. But given that I am using EPSG:32737 versus EPSG:4326, I want to know

A more precise definition of the problem is this. The original imagery is 6cm resolution, and that is in EPSG:32737. After converting the image to EPSG:4326, I want to resample the image to 1 meter resolution. So how many degrees correspond to 1 meter distance in Tanzania? Is there a simple equation to do this, or is there some package or website I can use?

In terms of tools, I was using rasterio and GDAL in python. I would like to do this through code instead of through a tool like QGIS.


texts: How to find ID's of nearest stations (points) in a specific radius

texts:
I am using ArcGIS 10.0 and I am doing rainfall data analysis.

I have a single layer including more than 1,000 stations (points) and I need name (STATION_NO) of neighbor stations in more than one specific distance for each station. I've tried all PROXIMITY options in Arctoolbox, but the outputs were not what I expected.

For example: I need to know the name of 5 nearest stations in a radius of 50Km.


texts: proj.4 definition for (obsolete) UK War Office "Italy South" grid

texts:
I'm helping a friend convert his late father's war diary coordinates to WGS84. I have a proj.4 definition that appears to work in QGIS, but if I transform a scanned map with gdalwarp, it's not even close. What am I doing wrong with my CRS definition, please?

Coordinates were given in the Modified British System, which is a variant of the OSGB grid, but extended over several zones in Europe. Projection parameters (given to me by the linked site's author) are as follows:

Projection - Lambert Conical Orthomorphic
Ellipsoid: Bessel 1841
False Easting : 700000
False Northing : 600000
Central Meridian : 14.0°
Central Parallel : 39.5°
Scale Factor : 0.99906

Everything except the scale factor can be confirmed from 100k Index to WWII topo maps of Italy at McMaster University. From this, I developed a proj.4 definition:

+proj=lcc +lat_1=39.5 +lon_0=14 +k_0=0.99906 +x_0=700000 +y_0=600000 +ellps=bessel +units=m +no_defs

After converting alphanumeric coordinates to grid coordinates and using the above definition, I get expected values for some sample locations:

rN175845 → 717500, 784500 → 41.1629 °N, 14.2086 °E
rN138862 → 713800, 786200 → 41.1783 °N, 14.1646 °E

When I try to use gdalwarp to translate a georeferenced map sheet to WGS84, however, I end up with latitudes around 3-4 °N. The command line I used is:

gdalwarp -s_srs "+proj=lcc +lat_1=39.5 +lon_0=14 +k_0=0.99906 +x_0=700000 +y_0=600000 +ellps=bessel  +units=m +no_defs" -t_srs EPSG:4326 -r cubic -of GTiff in_modified.tif out-4326.tif

gdalinfo reports the following for the input file:

Size is 6641, 5162
Coordinate System is `'
Origin = (694095.553313978714868,954540.625156613648869)
Pixel Size = (8.495191956798934,-8.495191956798934)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  694095.553,  954540.625)
Lower Left  (  694095.553,  910688.444)
Upper Right (  750512.123,  954540.625)
Lower Right (  750512.123,  910688.444)
Center      (  722303.838,  932614.535)
Band 1 Block=6641x1 Type=Byte, ColorInterp=Red
Band 2 Block=6641x1 Type=Byte, ColorInterp=Green
Band 3 Block=6641x1 Type=Byte, ColorInterp=Blue

(the sample coordinates are on another map sheet, btw.)

Update: a working proj.4 definition, based on mkennedy's comment below, is:

+proj=lcc +lat_0=39.5 +lat_1=39.5 +lon_0=14 +k_0=0.99906 +x_0=700000 +y_0=600000 +ellps=bessel +units=m +no_defs


texts: converting geojson to shapefile gives error "Unable to Unable to open datasource using following drivers"

texts:
Im trying to convert geojson file to shapefile. Below is the code

import ogr2ogr
import json

def main():
  #note: main is expecting sys.argv, where the first argument is the script name
  #so, the argument indices in the array need to be offset by 1
  ogr2ogr.main(["","-f","ESRI Shapefile","I:/shapefile/json_output.shp","I:/json_output.json -overwrite engine=GeoJSON"])

if __name__ == '__main__':
    main()

I get the following error

Unable to open datasource `I:/json_output.json -overwrite engine=GeoJSON' with the following drivers.
  ->  PCIDSK
  ->  netCDF
  ->  JP2ECW
  ->  JP2OpenJPEG
  ->  JPEG2000
  ->  DB2ODBC
  ->  ESRI Shapefile
  ->  MapInfo File
  ->  UK .NTF
  ->  OGR_SDTS
  ->  S57
  ->  DGN
  ->  OGR_VRT
  ->  REC
  ->  Memory
  ->  BNA
  ->  CSV
  ->  GML
  ->  GPX
  ->  LIBKML
  ->  KML
  ->  GeoJSON
  ->  OGR_GMT
  ->  GPKG
  ->  SQLite
  ->  ODBC
  ->  WAsP
  ->  PGeo
  ->  MSSQLSpatial
  ->  PostgreSQL
  ->  OpenFileGDB
  ->  XPlane
  ->  DXF
  ->  Geoconcept
  ->  GeoRSS
  ->  GPSTrackMaker
  ->  VFK
  ->  PGDUMP
  ->  OSM
  ->  GPSBabel
  ->  SUA
  ->  OpenAir
  ->  OGR_PDS
  ->  WFS
  ->  HTF
  ->  AeronavFAA
  ->  Geomedia
  ->  EDIGEO
  ->  GFT
  ->  SVG
  ->  CouchDB
  ->  Cloudant
  ->  Idrisi
  ->  ARCGEN
  ->  SEGUKOOA
  ->  SEGY
  ->  XLS
  ->  ODS
  ->  XLSX
  ->  ElasticSearch
  ->  Walk
  ->  Carto
  ->  AmigoCloud
  ->  SXF
  ->  Selafin
  ->  JML
  ->  PLSCENES
  ->  CSW
  ->  VDV
  ->  TIGER
  ->  AVCBin
  ->  AVCE00
  ->  HTTP

I tried the following link
Unable to open datasource 'MySQL' using ogr2ogr but no luck. Please help!!


texts: Has anyone any idea what the CRS is for France?  I am new to QGIS (moving over from using Arcgis) and trying to find out basics.  I have found a huge list beginning NTF (Paris) but wonder which one to apply.


texts: QgsMapToolEmitPoint - output coordinate to QGIS python console

texts:
I'm trying to output points to the QGIS Python Console with the following code which works:
canvas = iface.mapCanvas()
emit = QgsMapToolEmitPoint(canvas)
canvas.setMapTool(emit)

def display(pt, e):
    print(f"{pt} {e}")

emit.canvasClicked.connect(display)

However if I change
canvas.setMapTool(emit)

to
cc = QgsMapTool(canvas)
canvas.setMapTool(cc)

The maptool is set on the canvas but nothing prints to the console and I don't understand why.


texts: GetCoverage request - don't load result into memory and directly to hard drive

texts:
I'd be super grateful for any pointer to some resource, or any other kind of help. I want to make a GetCoverage-request against a server and download more or less 24 GB of data. In R I would do something like this:
library(httr)
res = GET("server...REQUEST=GetCoverage...")

My question now would be if I can download the data directly to my hard drive? Because this way, I'm assuming the data gets put into memory. At the moment the result is simply NULL and I don't know why. A smaller GetCapabilities and DescribeCoverage result in statuscode 200.


texts: Is it possible to show labels in QGIS CLOUD?

texts:
I have recently purchased the QGIS Cloud PRO version, but when I upload a map, I cannot visualize the labels of each layer.

Obviously, I have the labels activated on QGIS, but they just don't appear when uploading them to QGIS Cloud.
Is there any option in QGIS Cloud I am missing about labels?


texts: Definition Query Builder with multiple parts

texts:
I am working with agriculture and looking at two farms and using the query builder.

Each one has the same Farm Number (FSN) so I am trying to show particular tracts of land from those FSNs and particular Common Land Units (CLUs).

Is there a way to write a query to basically tell the system, "I want to see clus 1 and 2 from tract 443 and clus 1 and 4 from tract 442 both from farm 551?

I have been creating new layers for each tract, but I know there has to be a more elegant way to do this.

This is what I have for getting one tract:

state_code = '22' AND county_code = '053' AND farm_number = '516' AND tract_number = '228' AND clu_number IN ('17', '21')`

What would I need to do to get another tract and its clus from the same farm in this query?


texts: Output layer is not extent of both input layers in ArcMap

texts:
I have seen
Weighted overlay output raster does not match the extent of my input rasters
but it does not give me a solution.
I have two rasters and I'm using the raster calculator, the output it is giving me is the intersection, when I want the output to be spread across the full extent of both input layers (even when only one layer appears in that place).


texts: QGIS Georeferencer Plugin missing compression parameters

texts:
I am trying to georeference a rather large PDF. I was successful in doing so without compression.
However, when trying to make the map lighter, I added some compression. The options are: LWZ, PACKBITS, and DEFLATE. No matter which one I use, the output resolution is very low. The partial solution is checking Set target resolution, but still it only goes so far (1, -1), the output is pixeled and the labels on the old map I'm trying to georeference cannot be read.
I know this problem could be solved if I transform the pdfs into images and pre-resize them. But I'm looking for a more integrated solution. Is there a way to define the target resolution differently?
My setup:
QGIS: 3.10.10-A Coruña,
GDAL: 3.0.4, Win 10, FullHD laptop + 4K display.


texts:
I am trying to perform a classification in Arcmap's raster calculator, but it seems that the use of multiple conditional is confusing me. What am I doing wrong in the following con-statement?:

Con("raster" < 1,0.5, Con("raster" >= 1,1.5 & Con("raster" < 15,1.5, Con("raster" > 15,2.0))))


texts: Is there a vector layer showing sea level change by different cities around the world (so not just US or just EU ones)? I remember looking at ArcGIS Webmap some time ago, but can't find it now anymore. I would definitely prefer vector data as I would be able to have cities listed in the table, over raster data where I would need to manually look up for approximate number of each city.


texts: Stream Feature generated from DEM is missing data

texts:
I tried to search other questions to see if this has been answered but couldn't find anything, so sorry if this is redundant. I use the Hydrology toolset in Spatial Analyst to create a stream feature in ArcMap using a downloaded DEM. The series of tools I used was the following:

Fill (input DEM)
Flow Direction (input Fill)
Flow Accumulation (input Flow Direction)
Conditional > Con (input Flow Accum – input cond raster, Flow Dir input true rast & the optional)
Stream Order (input Con and Flow Dir)
Stream to feature (input Stream order & Flow Dir)
a.  Definition query: "GRID_CODE" > 2

I am relatively happy with the outcome except in the north section of the watershed where the stream vector detail was not read correctly and only produced straight lines.

Am I missing a step that would allow this to be read properly and converted into the same level of stream detail as everywhere else in the watershed?

Note: The straight lines are over cropfields, residences, and timber. I included the only waterbodies on the maps because I did notice on other stream generations that open water creates "straight lines." See aerial image below.


texts: How to Publish a Time Aware Mosaic dataset as an ImageService with Temporal component?

texts:
I have a Time Aware Mosaic dataset , which contains 5 rasters. I have entered the dates for the given rasters, by adding a Text Field with the dates. (The dates of the images are from 3/1/2013 to 3/5/2013)

The Mosaic dataset works great with the Time Slider window in ArcMap. But when I publish it to ArcGIS Server as an Image Service (through ArcCatlog), I do not see the correct Time Extent on the REST end point. The REST end point shows the Time Info as:

  Time Info:

  Time Extent: [2013/01/01 00:00:00 UTC, 2013/01/01 00:00:00 UTC]
     Time Reference: N/A

What steps have I missed? Where do you set the Time Extent/Time Reference while publishing the image service?

My Configuration: ArcGIS Dekstop 10.0 (ArcInfo)+ ArcGIS Server Advanced+ Image Extension


texts: Identifying slopes surrounding a lake (in ArcGIS)

texts:
I'm working with ArcGIS 10.7 to investigate lakes in the mountains. I have created a 2km "hazard buffer" around every lake-TIN showing instable slopes (>40°) above the lake surface from where a rockfall could hit the lake (see example below).

However, not all of the slopes inside the buffer are pointing towards my lake, some of them are on the other side of ridges, naturally.

I need to identify every mountainside with a slope >40° that is pointing towards my lake and classify them depending on their distance and height difference to the lake. The distance and height difference should not be too difficult to figure out (I hope) - but I struggle with the slope detection.

How could I use aspect or some other criteria to determine which slopes to keep?

EDIT: I know that my question is relatively broad. However, I'm not looking for a detailed workflow but for hints or tips on what I might look into.


texts:
Recently this community helped me to get a list of the tables from an MXD file:

IMapDocument mapDoc = new MapDocument();
mapDoc.Open("C://mymxd.mxd");

IActiveView act = mapDoc.ActiveView;
IMap map = act.FocusMap;

IStandaloneTableCollection coleccionTablas = (IStandaloneTableCollection)map;

for (int i = 0; i < coleccionTablas.StandaloneTableCount; i++)
{
    IStandaloneTable stTable = coleccionTablas.get_StandaloneTable(i);
    string nameTOC = stTable.Name;
}

That way I get the name of the table as shown in the "Table of Contents" of ArcMap. I thought that from there, it would be easy to get the name of the table as stored in the database. In other words: I need the datasource of the table.

I found this code: Changing data source of table using ArcObjects?

That same solution is the only one I have found so far in different sites. It seems pretty easy, but I only get ITable objects that are null... so getting the datasource for those is impossible.

This is my code:

for (int i = 0; i < coleccionTablas.StandaloneTableCount; i++)
{
    IStandaloneTable stTable = coleccionTablas.get_StandaloneTable(i);
    string nameTOC = stTable.Name;

    ITable table = stTable.Table; //--> table is always null!!
    IDataset dataset = (IDataset)table;
    IDatasetName datasetName = (IDatasetName)dataset.FullName;
}

I have tested several MXD files with tables, but the ITable object I get from the IStandaloneTable is always null and I have no idea why.

I have also tried with ITableCollection, but the result is the same.

Does anyone know why I keep getting nulls? Maybe there is something horribly wrong with my code.


texts: Automatically generate random polygons with a given number of vertices

texts:
Updated to make it clear.

What would be the fastest way to generate a random polygon in PostGIS with a specific number of vertices? I want to create a PostGIS function like

GeneratePolygon(1000, 4326)

That would generate a polygon of 1000 vertices using SRS 4326.

I'm not worried about the type of polygon, just about the number of vertices, for db performance. For example: will it be able to store and retrieve polygons with 100K vertices? How fast?

A possible way to generate polygons would be to sew 1000 points of a GPX track, forming a linestring, and make that the exterior bound of a polygon. I guess the function ST_MakePolygon could be used for something like that. The trick, I think, is how to generate the points. Something like the True part of this image


texts: Using Cell Size in raster calculator (ModelBuilder)?

texts:
I'm using ModelBuilder in ArcGis to create model which can calculate some stuff from my raster.

At some point I have to multiply raster value and raster cell size. I've created variable Cell Size, but it doesn't appear in Raster Calculator 'Layers and Variables'.

I could manually type value reffering to cell size, but it's not the point. In the future I'd like to use that model to make calculations on other raster with different cell size.

Does anyone know how to use Cell Size in raster calculator?


texts: Adding Icons to Configurable Map Viewer (CMV) Widgets?

texts: How do you add the Configurable Map Viewer (CMV) widget icons like the screenshot below?


texts: Seeking GeoJSON data for all counties of UK?

texts:
I'm having a look at http://code.highcharts.com/mapdata/countries/gb/gb-all.js

However I'm having trouble identifing a list of counties from this geojson data.

I've managed to get a rough list by checking for a property type feature of 'Administrative County'.


texts:
Is there a rubbersheet tool in QGIS similar to ArcMap?
Basically, I want the yellow line to move to the red line.

I have tried Vector Bender with a buffer of 0.1m through to 100m and still get a bad result:
Input

Bend Green to Red
Output


texts: Leaflet.geojson webpage performance out of hand, how to convert to .png?

texts:
I exported some countries as a single .geojson file. The file includes every country's vertices and coordinates, as that's what a .geojson is ─ a vector file.

Thing is, there are about 20 countries and considering each of them are 100-120 sided polygons, it's not surprising that the file size is +100 MB.

Leaflet has a geojson class that performs a geojson overlay for points, polygons etc. that reads the vector file and displays it on the map.

Since I want to post this .geojson file on the web using Leaflet, making every user download 100 MB for this is out of the question. Apart from that, when I tested this Leaflet map on my computer using Chrome, it takes about 5 seconds each time I zoom in/out for the map to come to resolution (the map crashes for 5 seconds because the webpage becomes unresponsive due to this vector file workload). I suppose pulling 100 MB from the RAM each time somebody zooms would be hard for just about any computer.

Is there any way to create a .png file of each country instead and overlay that onto the Leaflet map? I suppose that would greatly reduce file size and dramatically improve map performance.


texts: Selecting by location and attibutes simultaneously?

texts:
If I have a feature class ( say 'Cities') and I would like to select all features from another layer ('Incorporated Places') that nest within the cities layer and also have the same value for their 'Name' field.

How would I do this?

Specifically this will be used to remove duplicate labels between the incorporated place and its bounding city, and currently my method is to merge the centroids of the two features together and then remove duplicates between label classes, but I would prefer if there was something more elegant that would simply select the duplicate Place features instead.

I have seen other questions similar to this but the answers that were given for these suggest using ModelBuilder, which doesnt consider the fact that this process is based upon one selection with two conditions, rather than two separate selections, so as far as I  know ModelBuilder wouldn't be able to accomplish this.


texts: How do I style a GeoTIFF correctly in QGIS?

texts:
I downloaded the bathymetric data from GEBCO in geotiff format. When I load it in QGIS, the image is just completely grey. I am using QGIS 2.8.

How I can visualize this data in QGIS?


texts: "Create Random Points" does not accept integer column as "value from input field"

texts:
Using QGIS 2.6.1-Brighton in Ubuntu 14.04

I have a polygon layer loaded from a PostGIS query with int4 values that range from 0 to 3000. When I try using the Random Points tool on this layer, none of the integer columns appear in the dropdown list for "Use value from input field." When I divided the values in a column by 10 to diminish their range (while still using integers) I still could not run the function.

UPDATE:

I saved the query layer to file and reopened it and this still did not resolve the issue.


texts: Cannot retrieve attribute information from PostGIS layer

texts:
I want to develop a simple plugin and try to get started with Python in QGIS.

However, I have already run into troubles just loading my data properly.
From the cook book I am using the code below successfully to add my layer in QGIS. But it won't retrieve information from my layer:

uri = QgsDataSourceURI()
#set host name, port , database name, username and password
uri.setConnection("localhost", "5432", "local_db", "postgres", "postgres")
#set database schema, table name, geometry column and optionaly subset(WHERE clause)
uri.setDataSource ("friluftsliv", "climbing", "geom")
# Defining the layer name and layer type for QGIS?
climbing=QgsVectorLayer (uri .uri() ,"Climbing","postgres")
# Adding a layer to registry
QgsMapLayerRegistry.instance().addMapLayer(climbing)

Now I am trying to retrieve information from my layers using:

for field in climbing.pendingFields():
  print field.id_0()

This does not give me anything, no errors nor info... I have 27 rows in my table. It has SERIAL and PK in database.
 From http://www.qgistutorials.com/en/docs/getting_started_with_pyqgis.html
I tried with nothing to be printed.

for field in layer.getFeatures():
  print field['id_0']

I have another layer called hiking that I try to add with the same code as above. That layer don't get added to the registry. I don't get any errors when I defining the layer but if I run QgsMapLayerRegistry.instance().mapLayers() it does not come up in the list...

Any advice on what might be wrong?


texts: NetCDF to raster layer, spacing in their coordinate values

texts:
I'm trying to convert a NetCDF file to a raster layer in ArcGIS Pro, in order to make a space-time animation.
But I'm getting a error message that I don't know how to solve:
ERROR 000276: One or both dimensions have variable spacing in their coordinate values

I know that the spacing between x-coordinates must be equal and the spacing between y-coordinates must be equal, but the differences between latitudes decreases towards poles, so, how can I correct this?
This is the code that I have so far:
# Import system modules
import arcpy

# Set local variables
inNetCDFFile = "C:/Users/Miguel/Documents/ORIGINS/animationWRF/data/wrfout/wrfout_d03_2020-03-
09_00_00_00.nc"
variable = "CO2_ANT"
XDimension = "XLON"
YDimension = "YLAT"
outRasterLayer = "C:/Users/Miguel/Documents/ORIGINS/animationWRF/data/netcdf2raster/wrfout_d03_2020-
03-09_00_00_00"
bandDimmension = ""
dimensionValues = ""
valueSelectionMethod = ""

# Execute MakeNetCDFRasterLayer
arcpy.MakeNetCDFRasterLayer_md(inNetCDFFile, variable, XDimension, YDimension,
                           outRasterLayer, bandDimmension, dimensionValues,
                           valueSelectionMethod)

How can I solve this issue?


texts: How do I convert OSM XML to rendered tiles for offline use?

texts:
I'm new to OpenStreetMap and I'm trying to figure out how to process an OSM XML file and create map tiles for offline use. I already have the OSM XML file downloaded; I just don't know what to do next.

Once I have the tiles, I'll be displaying them in a control hosted in a Windows application. The basic application is complete. I'm now moving the data to local storage instead of online.

How do I take my OSM XML and process it into tiles?


texts: Adding Descriptions to custom ArcToolbox 9.3.1

texts: I looked http://webhelp.esri.com/arcgiSDEsktop/9.3/index.cfm?TopicName=Setting_script_tool_parameters and couldn't find how to add Descriptions for each parameter in the script tool.  I see how to add a description for the overall script tool but not for each parameter.


texts: Using Node.js to cluster millions of points on the server side

texts:
Situation:

I am testing out different visualization techniques for web mapping a data set that has about 2.5 million points and a lengthy set of attributes (almost 30 fields). So its quite a large file. I want to try clustering this data. I know sending that much data to a client will either result in a time-out or crash so I've been look at other methods and found the supercluster with mapbox. I used the top rated answer in this question to build a web map with Leaflet that handles about 500 000 points pretty well (approx 6 seconds to load). While this is good its only a subset of the data and its not as instantaneous as the Supercluster claims. Now I want to test using Node.js to do the clustering serverside, this way I can hopefully load the full 2.5 million point data set quickly.

Question(s):

I've never used Node.js before but I'm getting the gist of it. I need some help figuring out how to use the supercluster with Node.js because the GitHub site doesn't have clear instructions for beginners. I'm at the stage where I'm using the Demo index.html and its not showing me any points (I have not changed anything in the scripts).

Has anyone used the supercluster with Node before or has any general tips and tricks for me?

Also can I use GeoServer with Node.js?

Or if someone has a different suggestion for clustering 2.5 million points on the fly I'm interested in any and all options!

Specs:

OS -> Windows 10, Node.js -> 12.14.1, Leaflet -> 1.3.1


texts:
I built my own map canvas inside a plugin in QGIS and I want to add two layers (one by one) in order to overlay the two.
I have this code:

fileCI = "....\plugins\teste\test.shp"
fileInfo = QFileInfo(fileCI)
layer = QgsVectorLayer(fileCI, fileInfo.baseName(), "ogr")

if not layer.isValid():
    return

if self.canvas.layerCount()==0:

    # add layer to the registry
    QgsMapLayerRegistry.instance().addMapLayer(layer, False)

    # set extent to the extent of our layer
    self.canvas.setExtent(layer.extent())

    # set the map canvas layer set
    self.canvas.setLayerSet([QgsMapCanvasLayer(layer)])
    self.canvas.setVisible(True)

else:
    layers = self.canvas.layers()
    layers.insert(1,layer)

    # add layer to the registry
    QgsMapLayerRegistry.instance().addMapLayers(layers,False)

    for i in range(0,self.canvas.layerCount()+1):
        self.canvas.setLayerSet([QgsMapCanvasLayer(layers[i])])
        self.canvas.setVisible(True)

    for lay in layers:
        lay.triggerRepaint()

But when I add the first vector file, ok, when I add the second the map canvas stay empty.

How do I do to overlay the two layers?


texts: Rasterizing shapefile with rasterio and gdal gives all cells 0 values

texts:
I am trying to rasterize my shapefile (20 polygons) with rasterio and fiona and turn the new raster into a mask for clipping.

I need to rasterize my all polygons in my shapefile and turn on mask and clip my raster tif with mask  how can I do it ?

This is my first script, but that it give just value 0 all mask black.

import rasterio as rio
import fiona
import rasterio
import rasterio.mask

with fiona.open("C:/.../poly.shp", "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]

with rasterio.open("D:/.../MODIS_EVI_2006.tif") as src:
    out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
    out_meta = src.meta.copy()

out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})

with rasterio.open("mask2006.tif", "w", **out_meta) as dest:
    dest.write(out_image)

I'm also trying with gdal.RasterizeLayer and I had the same result.

from osgeo import gdal
from osgeo import ogr
from osgeo import gdalconst

ndsm = 'D:/../MODIS_EVI_2000.tif'
shp = 'C:/Doy2000/Pivos/1996poly.shp'
data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
#source_layer = data.GetLayer()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[3] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
mb_v = ogr.Open(shp)
mb_l = mb_v.GetLayer()
pixel_width = geo_transform[1]
output = 'C:/Doy2000/Pivos/my2000.tif'
target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))
band = target_ds.GetRasterBand(1)
NoData_value = -999999
band.SetNoDataValue(NoData_value)
band.FlushCache()
gdal.RasterizeLayer(target_ds, [1], mb_l, options=["ATTRIBUTE=hedgerow"])

target_ds = None