texts:
Is there any way to add multiple grids on the first parameter of saga:rastercalculator?

If I added one file as string at first parameters and the other files at second parameter as list, the algorithm worked, but it resampled my data and I don't need this.

I have some raster files with the same grid system generated in a for loop. I want  to sum their values in a single raster.

This is my code:

s1 = 'C:/Users/.../OUTPUT.sdat'
s2 = 'C:/Users/.../OUTPUT.sdat'
s3 = 'C:/Users/.../OUTPUT.sdat'
list = [s1,s2,s3]
rezult = processing.runalg('saga:rastercalculator', list, [], "g1+g2+g3", True, 8, None)

To simplify the code above, I used the files path as a string, but in my algorythm, the files paths are dynamically added to the list and the formula is dynamically updated.

When I run the above code, I get this error:

  Wrong parameter value: ['C:/Users/Flaviu/AppData/Local/Temp/processing6750b371ece74c718b442b56c8f3ef1f/c99bbf3a6dc64ba798285e21a5040648/OUTPUT11.sdat', 'C:/Users/Flaviu/AppData/Local/Temp/processing6750b371ece74c718b442b56c8f3ef1f/c99bbf3a6dc64ba798285e21a5040648/OUTPUT22.sdat', 'C:/Users/Flaviu/AppData/Local/Temp/processing6750b371ece74c718b442b56c8f3ef1f/c99bbf3a6dc64ba798285e21a5040648/OUTPUT33.sdat']

Is there an easier way to do this operation using QgsRasterCalculator ?


texts: How to install PostGIS correctly with OpenGeo Suite 4.0 community edition on Ubuntu?

texts:
I have installed the OpenGeo Suite community edition, version 3 and 4, onto Ubuntu Server 10.04 and 12.04 multiple times in the past, with virtually no problems.  With the latest version of OpenGeo Suite (which includes PostGIS 2.1), I am able to install almost all of the components successfully, except PostGIS, on ubuntu 10.04 and 12.04.  With earlier editions, PostGIS was installed as part of the full Opengeo installation process, but now, although PostGreSQL 9.3 is successfully installed, there are no PostGIS or opengeo database templates installed, nor are the spatial functions installed. According to apt-get, however, PostGIS 2.1 is installed.

I have also attempted to install OpenGeo v4 community edition on a local ubuntu server on a vm.  Same issue... geoserver, geowebcache, the opengeo dashboard, PostGres 9.3, and PostGIS 2.1 all installed, but no opengeo or postgis databases, and no functions.  I'm wondering if this is an issue with the Opengeo installation package?


texts: Wrong calculation of line start and end XY coordinates of a merged line- how can I fix this?

texts:
I have a feature class containing several rivers, which are derived from merging individual river segments to make a continuous line. They have been merged using the editor tool.
I have calculated the XY start and end coordinates for my merged lines but it is giving me random coordinates of positions which aren't the start and end of the line. Instead, it seems to be giving me start and end XY coordinates of pre-merged segments. How can I alter this so it recognises the line is now a continuous merged line and retrieve the XY start and end coordinates of the merged line not the individual segments.

See image attached. The coordinates are calculated as start and end XY using the calculate geometry tool but as you can see, they aren't at the start and end at all.


texts:
I have a raster and I'm trying to visualize it on QGIS.
For some reason I can only visualize the raster if I upload it alone without any other layer. If I add BING layer, it can't load due to python layer. If I upload first the BING and then the layer, the layer ''added'' but I can't see it.

I tried to save again my raster as WGS84 but  and also to register it with the BING and it didn't have any error, but still I can't see the results together, only if I upload it without the BING.


texts: Loading Satellite Imagery Into Trimble GeoExplorer

texts: How would I go about loading satellite imagery in TIFF format as a background file or image into my Trimble GeoExplorer?


texts:
I have a File Geodatabase Table(no geometry). I am trying to update the last column, called S_Name from concatenated values from other fields.
I am going to write out my logic so it is clear what I want to accomplish

if the PRE_TYPE field is not Null and the SUF_TYPE field is null then update S_Name to the values of PRE_TYPE+NAME+SUF_MOD

if PRE_TYPE is Null and SUF_TYPE is Not Null then update S_Name
to NAME+SUF_TYPE

Here is my code below. I am pretty sure the logic is sound. The script runs with no errors but does not print bah, and does not update the table.

import arcpy
arcpy.env.workspace = "C:path\\Tran_road_NJ.gdb"
layer = "C:path\\Tran_road_NJ.gdb\\Tran_road_geocode_NJ_1"
with arcpy.da.UpdateCursor(layer, ['PRE_TYPE','NAME','SUF_TYPE','SUF_MOD','S_Name']) as rd:
    for x in rd:
        if x[0] != None and x[1] != None and x[2] == None:
            bah = x[0]+' '+x[1] +' '+x[3]
            x[4] = str(bah)
            print bah
            rd.updateRow(x)
        elif x[0] == None and x[1] != None and x[2] != None:
            bah = x[1]+' '+x[2]
            x[4] = str(bah)
            print bah
            rd.updateRow(x)


texts:
I have a tiff raster with values from 0 to 255 (one band) and 7 colors. I would like to change colors of my tiff.

I tried with QGIS but I can't save it with new colors.

 from this

 to that

Somebody have an idea ?


texts: I can not download gpx from GPS through GPS tool with GPS babel on QGIS

texts:
I have already install GPS tool plug-in and GPS Babel.
Garmin GPS map 62s is connected my PC through USB cable.

I changed device setting according to below site.
http://nickmcw.wordpress.com/2010/07/30/quantum-gis-using-a-usb-gps-device/

I set below condition on "Download from GPS"tab.
GPS device: Garmin USB( ←I create by myself)
Port:usb
Feature type:way point
Layer name:F:\Garmin\GPX
output file:test

At least, I got new GPX file on layer of QGIS.

Altgough the new GPX file is expected include existing gpx information from GPS, there is no row in attribute table.

Colum names were put automatically.followings are.
Name/elevation/symbol/comment/description/source・・・・・

Why I can not get gpx file from GPS???


texts: GRASS session fail to start because of g.extension failed on linux mint

texts:
I am struggling to start any existing Grass sessions (that was working before) in Linux Mint 18.3 Cinnamon 64 bit on an Oracle virtual machine.

The following error is popping up every time I am trying to open any previously working GRASS session.

I have tried the following:

Uninstall GRASS 7.2 and all other grass components (GRASS Core, GUI)
and reinstalled it.
Restarted the Linux after uninstalling GRASS
Double checked that no other GRASS session is open, or existing lock
is not there

Still cannot use GRASS. Any idea??


texts:
In my database, I've found a table (SDE_GEOMETRY276) with permissions granted to one user. My first thought was the "276" was either a layer_id or table registration_id. Searching sde_layers and sde_table_registry shows nothing with an id of 276.

Is it possible this is an orphaned table and, if so, what's the safest way of getting rid of it?


texts: How do I deal with tiles in the Orfeo Toolbox segmentation?

texts:
I'm using the Orfeo Toolbox for segementation of large Very High Resolution Satellite images. It performs very well and i'm mostly very happy with it. But I still get some weird/bad results that I assume are related to the partition of the process in tiles.

I'm doing this:

#Segmentation variables
    spatialr = '7'
    maxiter = '6000'
    ranger = '2,6'
    minsize = '8000'
    tilesize = '1048'
    simplify = '0.4'
    # SEGMENTATION
    print "Performing image segmentation"
    app = otbApplication.Registry.CreateApplication("Segmentation")
    app.SetParameterString("in", imagem_base_com_ndvi)
    app.SetParameterString("filter.meanshift.spatialr", spatialr)
    app.SetParameterString("filter.meanshift.maxiter", maxiter)
    app.SetParameterString("filter.meanshift.ranger", ranger)
    app.SetParameterString("filter.meanshift.minsize", minsize)
    app.SetParameterString("mode.vector.inmask", image_mask)
    app.SetParameterString("mode.vector.stitch", 'True')
    app.SetParameterString("mode.vector.out", vector_out)
    app.SetParameterString("mode.vector.simplify", simplify)
    app.SetParameterString("mode.vector.outmode", 'ovw')
    app.SetParameterString("mode.vector.tilesize", tilesize)
    app.SetParameterString("mode.vector.neighbor", 'True')
    app.ExecuteAndWriteOutput()

Here are some examples of the resulting segmentation:

Anyone has ideas to help me avoid these errors?

Thanks.


texts: Looking for map calibration software and city database

texts:
I'm in need of a mapping software which, given a map, lets you calibrate it (choose which geographic coordinates project to which pixels on the map image) and then given a database of cities, you can navigate to the specified city on the map you just calibrated.

To reword the question, I need a mapping software that can calibrate my map image and accepts some kind of database format to display cities on that map.

Where to find such compatible map software and database? Preferably free and open source software and data.

If the question is too broad, I appreciate any points making me near to the answer.


texts: ArcGIS Route (Change direction phrase structure)

texts: Is it possible to change direction text phrasing? at the moment on the route directions, it mostly starts with Turn(right/left) depending on the esriDMT maneuverType. But some people prefer to see an instruction started with "continue then turn (right/left)" rather than "turn(right/left) then continue"?


texts: ArcGIS Server 'HTTP and HTTPS' security switching back to HTTP on reboot

texts:
I have an Amazon EC2 instance running ArcGIS server 10.3.1. It has been running with security set https. I recently had to change security to mixed mode http and https. Now when the server reboots security changes back to http only, breaking the web adapter that is configured to use SSL. After each reboot I have to change security back to mixed.

Before the change to mixed security the server was running on https only without problem for several years.

Is this expected behavior and is there a way to force the instance to stay on mixed security?


texts: Iterate over all attributes in QGIS Graphical Modeler

texts:
I would like to iterate over N unknown attributes in the QGIS Graphical Modeler.
For example dividing all attributes in an input layer by its area to create an area weighting.

I realise it could be done with Python, but I'm wondering if there is some clever way to accomplish this with out-of-the-box tools.


texts: arcpy.mapping Graduated Colors Symbology

texts:
I'm trying to figure out how to apply graduated color symbology to a feature class by using arcpy.mapping in ArcMap 10.1.

To start, a script runs and a feature class is added to a map document, no problem. This feature class needs the graduated colors. The value field used for the graduated colors is created during the script.

In order to apply the graduated colors, I have to create a layer file from the feature class. This is where I start to get confused. I know UpdateLayer has to be used, and I have to base the update on a source layer? I tried this with no luck.

If someone can just guide me on what order the steps are to apply a graduated color symbology to a feature class using arcpy.mapping, it will be appreciated. I should be able to figure out the rest of the scripting. Thanks!


texts: What is map file in GRASS GIS addon r.fuzzy.system and how it can be accessed?

texts:
I am inexperienced with GRASS GIS and I am trying to reproduce the fuzzy systems manual example. I am failing to load the MAP file, even if I wrote the 'txt' file with the rules and maps name.
Here is the manual example:https://grass.osgeo.org/grass78/manuals/addons/r.fuzzy.system.html.
What is a 'MAP file', how can it be accessed, and how do I use it?


texts: Should states/provinces be stored as geographic or geometric data types?

texts:
While reading the GeoDjango documentation, they advise using geometric rather than the geographic data type when storing values that are on a continental scale.

The Canadian province of Ontario, for example, is fairly large. Although it's not on a continental scale, will I face a significant loss of accuracy when performing polygon calculations and mapping using geometric (as opposed to geographic) storage? At what threshold should one switch over from using geometric data types to using geographic data types?


texts: Exporting basemap tiles Contextily (ctx.bounds2raster) gives a black/empty GEOTIFF

texts:
When exporting a tile basemap (osm topo tile service) from Contextily using ctx.bounds2raster I get a black geotiff raster with no data in it. The properties of the raster say it's an 8 bit unsigned integer- which is surely not compatible with a 3 band raster representing the tile images.
The raster seems to be missing the CRS as well, but when I set the layer CRS to 3857, it projects to the correct place at least.
There is a CRS related error:
'PROJ: proj_create_from_database: Cannot find proj.db
-87658.19839440205 -65044.452373738226 7082306.226668285 7109249.279867188
Error  [line 29] The WKT could not be parsed. OGR Error code 6
Here's the code:
import matplotlib.pyplot as plt
import geopandas as gpd
import geopy
import contextily as ctx
from shapely.geometry import Point, Polygon, LineString
import pandas as pd
import sys, traceback

def excep():
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]
    errorM =  str(sys.exc_info()[1]).split("\n")[0]
    if "," in errorM:
        errorM = errorM.split(",")[0] + ". " + errorM.split(",")[1]
    lineN =  str(tbinfo).split(", ")[1]
    mess = "Error  [" + lineN + "] " + errorM
    print(mess)

def add_basemap(ax, zoom, source):
    try:
        xmin, xmax, ymin, ymax = ax.axis()
        print (xmin, xmax, ymin, ymax)
        basemap, extent = ctx.bounds2img(xmin, ymin, xmax, ymax,  ll=False, zoom=zoom, source=source)
        ax.imshow(basemap, extent=extent, interpolation='bilinear')
        basemap, ext = ctx.bounds2raster(xmin, ymin, xmax, ymax, r"C:\Temp\topo.tif",zoom=zoom)
    except:
        excep()

fig, ax = plt.subplots(1, figsize=(8, 5))

minx = 481062
maxx = 493021
miny = 405741
maxy = 420049
epsg = 27700

# Create a simple vector bounding box for my area and plot with geopandas
polyList = [(minx,miny),(minx,maxy),(maxx,maxy),(maxx,miny),(minx,miny)]
boundingbox_geom = Polygon(polyList)

area_df = pd.DataFrame(data= [['box',boundingbox_geom]],columns=['box','geometry'])
area_gdf = gpd.GeoDataFrame(area_df,crs='epsg:{}'.format(epsg),geometry='geometry')

# get vector data into same CRS as tile services - Web Mercator
area_gdf2 = area_gdf.to_crs(epsg=3857)
area_gdf2.plot(ax=ax,alpha=0.5,color='yellow',lw=2,zorder=2,label='box')

topo = "https://a.tile.opentopomap.org/{z}/{x}/{y}.png"
add_basemap(ax, zoom=12,source=topo)
plt.show()

It plots ok in web mercator (you can ignore the yellow box as that's just me adding some vector data using geopandas):

But the geotiff is black/empty (the black box matches the tiled basemap in the plot above):

May be there are some syntax errors or I've misunderstood the contextily methods?
I would ideally like to request the tiles in EPSG:27700 so I can keep my plot in the CRS I want it in but wasn't sure how to do that, something like source_url&epsg=27700... but I will definitely need to export the image to a geotiff anyway for other reasons.


texts:
I just upgraded from QGIS 2.18 to 3 on my MacBook Pro (retina display, running High Sierra).  Some of the cursors are enormous and I can't figure out how to change their size.  For example, the Pan Map cursor is the normal size, but the Zoom In/Zoom Out cursors are about 2 or 3 times the normal size.  Screen shots show what I mean.

I don't have this problem with my iMac, so I'm wondering what's wrong on my MacBook Pro.  I don't think it's a retina screen issue, because it didn't have this issue with QGIS 2.18.

Any ideas?


texts: Row ID Column Type : User Maintained - Insert using IFeatureBuffer - ArcGIS 10.2

texts:
I am trying to insert a record in SDE feature class (which is not registered with SDE) using Arcobjects 10.2 using IFeatureBuffer.InsertFeature functionality. This throws an error message "The current table doesnt have an ArcSDE-maintained rowid column.".

The same code snippet is working fine using Arcobjects 10.0.

Questions:

a) Has IFeatureBuffer interface has changed? ( I don't see any updates in help documentation in ESRI )
b) Is there new rule imposed in SDE to insert a data like Row ID column type should be SDE maintained in 10.2?
c) Is there any impact on altering existing feature class to change the rowid type?
d) Command to alter the feature class rowid column type property?

Attribute Administration Utility

Created by layer gen

Table Owner : PROJECT_SCHEMA
Table Name : TABLE1_
Registration Id : 7
Row ID Column : ID
Row ID Column Type : User Maintained
Row Lock :
Minimum Row ID :
Dependent Objects : Layer
Registration Date : Sun Sep 25 20:13:28 2011
Config. Keyword : DEFAULTS
User Privileges : SELECT, UPDATE, INSERT, DELETE
Visibility : Visible


texts: Georeferencing a CAD file and converting to SHP in ArcGIS without distance distortion

texts:
i have a land parcels dwg file in an unknown co-ordinate system and i want to georeference the file in ArcMAp and convert it to shapefile but the problem is that whenever i georef the cad drawings directly or or through spatial adjustment by converting it into shapefiles, the distances get distorted. for instance, if the length of line is 31 metres it changes to 0.000334 metres after projecting and georeferencing drawing.  i have my control points in UTM easting and northing. can anyone tell me the detailed method to finish this task?
The precise steps that i performed are:

importing dwg file in arcmap and assigning projected co-ordinate system to file.
georeferencing the dwg file through two control points by using an already georefenced tif file for destination points.

my cad drawing overlays on the tif file but main problem are the distances an area. they are inaccurate.


texts:
I'm running a ArcGIS 9.3 Oracle SDE geodatabase.
Our ArcMap client is extremely slow at displaying versionned feature class.
One table seems to be impacted : sde.state_lineages
A simple SELECT * FROM sde.state_lineages takes more than 10 seconds to display 69 rows.
Any idea what could be the problem?


texts:
Is there any freely available pointcloud data?

I'm especially interested in dense clouds of small areas(5m to 100m radius) but anything else would be fine too.

I've found that these pages have very nice and easy to access point cloud data.

http://kos.informatik.uni-osnabrueck.de/3Dscans/
http://www.liblas.org/samples/
http://www-graphics.stanford.edu/data/3Dscanrep/


texts:
My office will be seeing a big change in its GIS section. This section has been operational since the 1980's and has a huge collection of GIS data (i.e., shapefiles, raster files, data, etc) but never been through any inventory. Now it will happen.

Is there any automated way to extract all the information about the GIS data (i.e., shapefile, arc-info coverage, layer file, *.mxd, gdb, raster file, and more) from a PC to an Excel file? The information may include date of creation, date last edited, folder or container name, etc.


texts: Topographic correction (for land use classification) with SAGA or GRASS

texts: Is there a module in GRASS or SAGA (or a alternative open program) to do a topographic correction of a landsat image to enhance the land use classification?


texts: Calculate best area with overlaying multiple layers buffer shapes

texts:
I am using QGIS and I have a set of 20 layers with buffers. I'm interested in where they overlap most and I want to create a single layer that shows which are the places where they overlap most and how much. Somehow I already achieved the visual result changing the color and the transparency of the layers but I want just one layer that tells me how many layers are overlapping on a specific area.

Is it possible?


texts: django.contrib.gis.gdal.CoordTransform does not work

texts:
I want to convert EPSG:7844 toEPSG:4326 by using django.contrib.gis.gdal. I define CoordTransform in this way:
from django.contrib.gis.gdal import CoordTransform, SpatialReference
ct = CoordTransform(SpatialReference(7844), SpatialReference(4326))
print(ct)

Output:
Transform from "GDA2020" to "WGS 84"

Then I define point
from django.contrib.gis.geos import Point
point = Point(335810.4049635191, 6235191.808044555, srid=7844)
print(point)

Output:
SRID=7844;POINT (335810.4049635191 6235191.808044555)

And when I apply the transformation, it seems it does not work:
point.transform(ct)
print(point)

Output:
SRID=4326;POINT (335810.4049635191 6235191.808044555)

Does django.contrib.gis.gdal support such transformations? If yes, what am I missing here?


texts:
In ArcCatalog, getting the properties on an Annotation Feature Class exposes an Annotation Classes tab and within that is a list of the names of the associated annotation classes.

Assume an MXD with a feature-linked annotation layer and another layer with the associated source feature class, both stored in SDE. Can someone provide some code that will get that list of annotation class names?   I've been searching through the documentation with no success.  I'm using C#.


texts: QGIS extracted azimuth used for symbol rotation seems not reproduced faithfully

texts:
A field of calculated azimuths* (degrees) of digitized two-point lines was created. Field expression to calculate azimuth:

degrees( azimuth( point_n( $geometry ,1) , point_n( $geometry ,2 ) ) ) .

Using Vector Geometry / Extract Specific Vertices using node 1, a point layer was created from the line layer. The point layer included the azimuth field and the process created an angle field with identical values as the calculated azimuths. For the point layer, a linear-based svg symbol using a data-defined override of the azimuth for rotation (degrees required), symbols appear close to original only when near azimuths of 90 and 180 degrees but are off by what appears to be up to 15 or 20 degrees at intermediate rotations. Has anyone encountered and resolved this?


texts: Seeking Metropolitan Statistical Area (MSA) KML data?

texts:
GIS is not my field. I need this data for an odd job at work.
I have been unsuccessful in finding the data from Google FusionTables but found (what I believe to be) the shapefile for MSA's from the US Census on the US Census site
Is there an easy way to convert this file to KML polygon data?
I have tried QGIS but it seems to throw an error when I save the Layer as a KML file type.
Alternatively, is there any place I can simply find this polygon data?
Here is the error text I receive from the QGIS:
"Export to vector file failed.
Error: creation of data source failed (OGR error:Failed to create KML file polygon_data.kml.)"

@mtn.biker successfully created the file for me and outlined their process in the comment below if anyone else runs into this issue.


texts: How to test polar stereographic projection for a specific longitude intersection?

texts:
I'm writing a library that tests grids defined by map projections for intersections, and am running into an issue implementing the test for polar stereographic and rectangular projections - 'rectangular' meaning a grid where the sides of the projected grid coincide with latitude/longitude lines (e.g. mercator).

Without getting into the math too much, the formula to convert a lat/long coordinate to a north polar stereographic grid coordinate is:

E = (1 + sin(DxDyLatitude)) * EarthRadius
R = E * cos(Latitude) / (1 + sin(Latitude))
X = (sin(Longitude - Orientation) * -R) / IncrementX
Y = (cos(Longitude - Orientation) * R) / IncrementY

IncrementX/Y** - the distance in meters between grid units of X/Y
DxDyLatitude** - the latitude at which IncrementX/Y are valid
Orientation** - the longitude that coincides with the +Y axis of the stereographic projection
X/Y - the grid coordinate of the lat/long point

(** constant for a given grid)

To determine whether the top or bottom latitude of the rectangular grid intersects the stereographic grid, fix either X or Y at a grid boundary (left/right or bottom/top) and solve for the other.  If the other is within left/right or bottom/top and the longitude at the same point is within the rectangular grid, the grids can be said to intersect.  This is done by re-arranging the above equations:

Longitude = asin((X * IncrementX) / -R) + Orientation
Longitude = acos((Y * IncrementY) / R) + Orientation

However, determining whether the left or right longitude of the rectangular grid intersect seems not as easy due to the equation for R.  When re-arranging the equations to solve for latitude instead of longitude, I am stuck at:

cos(Latitude) / (1 + sin(Latitude)) = X / ((sin(Longitude - Orientation) * -E) / IncrementX)
cos(Latitude) / (1 + sin(Latitude)) = Y / ((cos(Longitude - Orientation) * E) / IncrementY)

I don't know of any trigonometric identity that reduces cos(x) / (1 + sin(x)), or if the equation can be re-arranged so that there is such an identity, or of any other mathematical equation for stereographic projections that will allow me to do this.


texts: Fixing document path during PDF export from ArcMap?

texts:
The document path on the side of my map changes when I export to PDF.
The text says: "dyn type="document" property="path"
It displays the filepath correctly when I'm editing the document (maps are located on a shared server).
However, when I export it to a PDF, in the PDF it shows the file as being in a temp folder on my local hard drive, a folder which doesn't even exist.

Why does the filepath switch when exported and how do I fix this?


texts:
I'm adding a layer with time control throught the widget "add Data"

But this layer doesn't show in map. Could someone give me a clue as to why this happens?


texts: How can I conduct dimension reduction using Geographically Weighted Principal Component Analysis?

texts:
I want to know how to compute the importance of components for geographically-weighted-principal-component-analysis(gwpca)? That is, I want to reduce dimension according to cumulative proportion of component.

Here is the code about pca and gwpca from the jstatsoft paper about GWmodel package

require(GWmodel)

data("DubVoter")

Data.scaled <- scale(as.matrix(Dub.voter@data[,4:11]))

pca.basic <- princomp(Data.scaled, cor = FALSE)
(pca.basic$sdev^2 / sum(pca.basic$sdev^2))*100

pca.basic$loadings

R.COV <- covMcd(Data.scaled, cor = FALSE, alpha = 0.75)
pca.robust <- princomp(Data.scaled, covmat = R.COV, cor = FALSE)
pca.robust$sdev^2 / sum(pca.robust$sdev^2)

pca.robust$loadings

Coords <- as.matrix(cbind(Dub.voter$X, Dub.voter$Y))
Data.scaled.spdf <- SpatialPointsDataFrame(Coords, as.data.frame(Data.scaled))

bw.gwpca.basic <- bw.gwpca(Data.scaled.spdf,
                           vars = colnames(Data.scaled.spdf@data), k = 3, robust = FALSE, adaptive = TRUE)
bw.gwpca.basic

bw.gwpca.robust <- bw.gwpca(Data.scaled.spdf,
                            vars = colnames(Data.scaled.spdf@data), k = 3, robust = TRUE, adaptive = TRUE)
bw.gwpca.robust

gwpca.basic <- gwpca(Data.scaled.spdf,
                     vars = colnames(Data.scaled.spdf@data), bw = bw.gwpca.basic, k = 8,
                     robust = FALSE, adaptive = TRUE)

gwpca.robust <- gwpca(Data.scaled.spdf,
                      vars = colnames(Data.scaled.spdf@data), bw = bw.gwpca.robust, k = 8,
                      robust = TRUE, adaptive = TRUE)

prop.var <- function(gwpca.obj, n.components) {
    return((rowSums(gwpca.obj$var[, 1:n.components]) /
                rowSums(gwpca.obj$var)) * 100)
}

var.gwpca.basic <- prop.var(gwpca.basic, 3)
var.gwpca.robust <- prop.var(gwpca.robust, 3)

Dub.voter$var.gwpca.basic <- var.gwpca.basic
Dub.voter$var.gwpca.robust <- var.gwpca.robust

mypalette.4 <-brewer.pal(8, "YlGnBu")

spplot(Dub.voter, "var.gwpca.basic", key.space = "right",
       col.regions = mypalette.4, cuts = 7,
       main = "PTV for local components 1 to 3 (basic GW PCA)",
       sp.layout = map.layout)

spplot(Dub.voter, "var.gwpca.robust", key.space = "right",
       col.regions = mypalette.4, cuts = 7,
       main = "PTV for local components 1 to 3 (robust GW PCA)",
       sp.layout = map.layout)

loadings.pc1.basic <- gwpca.basic$loadings[,,1]
win.item.basic <- max.col(abs(loadings.pc1.basic))

loadings.pc1.robust <- gwpca.robust$loadings[,,1]
win.item.robust <- max.col(abs(loadings.pc1.robust))

Dub.voter$win.item.basic <- win.item.basic
Dub.voter$win.item.robust <- win.item.robust

mypalette.5 <- c("lightpink", "blue", "grey", "purple",
                 "orange", "green", "brown", "yellow")

spplot(Dub.voter, "win.item.basic", key.space = "right",
       col.regions = mypalette.5, at = c(1, 2, 3, 4, 5, 6, 7, 8, 9),
       main = "Winning variable: highest abs. loading on local Comp.1 (basic)",
       colorkey = FALSE, sp.layout = map.layout)

spplot(Dub.voter, "win.item.robust", key.space = "right",
       col.regions = mypalette.5, at = c(1, 2, 3, 4, 5, 6, 7, 8, 9),
       main = "Winning variable: highest abs. loading on local Comp.1 (robust)",
       colorkey = FALSE, sp.layout = map.layout)

The pca.basic indicates we can keep the first to 8th component because of cumulative proportion of comp.8  is higher than 85%. But how can I conduct similar job with gwpca ?

Some related materials :

How can I conduct Geographically Weighted Principal Component Analysis using ArcGIS, Python and SPSS/R?

Harris P, Clarke A, Juggins S, et al. Enhancements to a Geographically Weighted Principal Component Analysis in the Context of an Application to an Environmental Data Set[J]. Geographical Analysis, 2015, 47(2): 146-172.


texts: ArcObjects Maplex labeling doesn't work?

texts:
I am trying to modify a FeatureLayer's Maplex label settings. For some reason those properties don't get modified when opening the document afterwards.
The document I am opening HAS already set the option "Use maplex label engine" property set.

def ApplyModificationsMainLayer(featureLayer):
    featureLayer.Name = "PAD"
    pGeoFeatureLayer = CType(featureLayer, esriCarto.IGeoFeatureLayer)
    pGeoFeatureLayer.DisplayAnnotation = True
    propertiesCollection = pGeoFeatureLayer.AnnotationProperties

    pLabelEngineLayerProperties = NewObj(esriCarto.MaplexLabelEngineLayerProperties, esriCarto.ILabelEngineLayerProperties2)

    pMaplexOverposter = NewObj(esriCarto.MaplexOverposterLayerProperties, esriCarto.IMaplexOverposterLayerProperties)
    pMaplexOverposter.FeatureType = 1
    pMaplexOverposter.LinePlacementMethod = 7
    pMaplexOverposter.PrimaryOffset = 10
    pMaplexOverposter.PrimaryOffsetUnit = 0
    pMaplexOverposter.CanReduceFontSize = True
    pOverposter = CType(pMaplexOverposter, esriCarto.IOverposterLayerProperties)
    pLabelEngineLayerProperties.OverposterLayerProperties = pOverposter

    pLabelEngineLayerProperties.Expression = "[Length]"
    pLabelEngineLayerProperties.Symbol.Size = 7

    propertiesCollection.Clear()
    propertiesCollection.Add(pLabelEngineLayerProperties)
    #cast to ILayer before saving
    pLayer = CType(pGeoFeatureLayer, esriCarto.ILayer)
    return pLayer

I have tried to get a reference of the "Default" Label Class and it seems that when I try to cast it's OverposterLayerProperties to IMaplexOverposterLayerProperties it returns a Null Type Pointer.

Why would this happen ?


texts: Difference in area of polygon and excised raster layer

texts:
I'm using QGIS 3.8. I have dowloaded GLOBCOVER 2009 raster land cover map and ecoregions shapefile. Polygon of ecoregions are used for mask to cut raster file. On new layer r.report was performed to measure land cover types.

I have found difference between area of polygon and excised raster layer. For example, Ural montane forest and tundra ecoregion is 175533 km2 and the result of r.report:

Total area is 116576,1 km2. Difference is 58956,88 km2. What causes the difference?


texts: Getting index of virtual field using PyQGIS

texts:
Im using QGIS 3.18.2 and I want to get the expression string of a virtual field in PyQGIS.
QgsVectorLayer.expressionField() is supposed to return the expression string. It takes the index (index: int) of the expression (virtual) field as argument.
Let's say I have a virtual field named expr_f_1. When I call  QgsVectorDataProvider.fieldNameIndex("expr_f_1") it returns -1, which corresponds to a missing field. This seems to happen for virtual fields in general.
They also don't appear in QgsVectorDataProvider.fieldNameMap().
How do I get the correct index of a virtual field, in order to pass it to  QgsVectorLayer.expressionField() to return the corresponding expression string?


texts: Parcel Fabric vs. Polygon feature class (for land parcel management)

texts: I'm wondering if I should migrate my existing parcel polygons (stored as a feature class in a geodatabase) to a parcel fabric in ArcMap. What are the advantages of doing this? I work for a county so the parcels are used mainly to store attributes (who owns which parcels in a graphic representation etc.) and general information - NOT for being survey accurate.


texts:
How can I access osm2po Webservice in my local network?

I can access http://localhost:8888/Osm2poService but if I try (on the same machine) to use the assigned ip address with http://192.168.1.5:8888/Osm2poService it doesn't work. Is there any configuration I need to change?


texts: Geoserver feature extraction "Failed to retrieve value for input features"

texts: I am trying to run the process gs:BufferFeatureCollection in geoserver Demo builder but i am getting the error "gs:BufferFeatureCollectionBufferBuffers features by a distance value supplied either as a parameter or by a feature attribute. Calculates buffers based on Cartesian distances.Failed to retrieve value for input features Unsupported class file major version 57"


texts: ArcGIS 10.1 Legend - Increase Icon on legend without increase on map

texts:
I was wondering if anyone knows how to increase an icon size on the legend without it increasing in size on the map. I have to print at A0, and when I look at the legend the icons for each type are quite small. Any attempt at increasing the size also increases it on the map.

I have tried what MAJ742 suggested - the legend does increase in size - but not the icons see attached image. I was looking to be able to enlarge the icons only on legend.


texts: Using Gdalwarp to Change Resolution of File in Degrees to Metres

texts:
I am using the command: gdalwarp -t_srs WGS84 -te 144 -36 150 -32 -crop_to_cutline AustTopo250_SI55.tif

I need to change the pixel resolution to 50m x 50m but the option is given as
 -tr xres yres: set output file resolution (in target georeferenced units).

The target unit is degrees. So the question is: How do I generate the output file as 50m x 50m pixels?

Feel free to explain relationship between metres and degrees.

Target File:

E:\>gdalinfo AustTopo250_SI55.tif
Driver: GTiff/GeoTIFF
Files: AustTopo250_SI55.tif
Size is 18750, 12500
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (144.000000000000000,-32.000000000000000)
Pixel Size = (0.000320000000000,-0.000320000000000)
Metadata:
  AREA_OR_POINT=Area
  COLORSPACE=RGB
  COMPRESSION_RATE_TARGET=9
  VERSION=2
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 144.0000000, -32.0000000) (144d 0' 0.00"E, 32d 0' 0.00"S)
Lower Left  ( 144.0000000, -36.0000000) (144d 0' 0.00"E, 36d 0' 0.00"S)
Upper Right ( 150.0000000, -32.0000000) (150d 0' 0.00"E, 32d 0' 0.00"S)
Lower Right ( 150.0000000, -36.0000000) (150d 0' 0.00"E, 36d 0' 0.00"S)
Center      ( 147.0000000, -34.0000000) (147d 0' 0.00"E, 34d 0' 0.00"S)
Band 1 Block=18750x1 Type=Byte, ColorInterp=Red
  Description = Red
Band 2 Block=18750x1 Type=Byte, ColorInterp=Green
  Description = Green
Band 3 Block=18750x1 Type=Byte, ColorInterp=Blue
  Description = Blue

Source File:

E:\>gdalinfo 76943_Aust_gda94.ecw
Driver: ECW/ERDAS Compressed Wavelets (SDK 3.x)
Files: 76943_Aust_gda94.ecw
Size is 128490, 109375
Coordinate System is:
GEOGCS["GEOCENTRIC DATUM of AUSTRALIA",
    DATUM["GDA94",
        SPHEROID["GRS80",6378137,298.257222101]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (112.883360000000040,-8.999920000000039)
Pixel Size = (0.000320000000000,-0.000320000000000)
Metadata:
  COLORSPACE=RGB
  COMPRESSION_RATE_TARGET=9
  VERSION=2
Corner Coordinates:
Upper Left  ( 112.8833600,  -8.9999200) (112d53' 0.10"E,  8d59'59.71"S)
Lower Left  ( 112.8833600, -43.9999200) (112d53' 0.10"E, 43d59'59.71"S)
Upper Right ( 154.0001600,  -8.9999200) (154d 0' 0.58"E,  8d59'59.71"S)
Lower Right ( 154.0001600, -43.9999200) (154d 0' 0.58"E, 43d59'59.71"S)
Center      ( 133.4417600, -26.4999200) (133d26'30.34"E, 26d29'59.71"S)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Description = Red
  Overviews: 64245x54687, 32122x27343, 16061x13671, 8030x6835, 4015x3417, 2007x1708, 1003x854, 501x427, 250x213
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Description = Green
  Overviews: 64245x54687, 32122x27343, 16061x13671, 8030x6835, 4015x3417, 2007x1708, 1003x854, 501x427, 250x213
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Description = Blue
  Overviews: 64245x54687, 32122x27343, 16061x13671, 8030x6835, 4015x3417, 2007x1708, 1003x854, 501x427, 250x213


texts: How to edit attribute tables using Python/ArcPy?

texts: I am new to Python, so sorry for such a trivial question. There is a shapefile (let's call it Structures.shp) and the script I'm writing needs to update the value of only one record in only one field (let's call the field structuretype). The field is the subtype field and the default value is (for example) 3, but the script needs to change the value for this particular record to 4. What is the easiest way for the script to update the value of this one cell?


texts: Finding number of lines passing through common segment of lines from overlapping multiple line shapefiles?

texts:
I have 16 different origin-destination route’s shapefile. 16 line shapefiles are different. I want to overlay these lines and find the common segments of these overlays routes. Suppose, one line is A-B, another is P-Q, another one is X-Y; I want to know how many routes go through C-D segment, how many go through E-F segment. I mean I want to separate the number of the routes by link.

How can I do this?

I imported the shapefiles from google map and shapefile have no feature class.


texts:
I have two points in OSGB 1936 (EPSG 27700):

A: 308188.48 608846.16

B: 308288.48 608846.16

Distance Matrix produces a result between A and B which is not exactly 100m (100.03935792736038 to be precise!). Why is that?


texts: Change boolean value in QGIS by CASE WHEN

texts:
I'm trying to change a field value by the next conditional SQL script from a PostGIS database in QGIS Field Calculator:

case
when "Field_1" = 'YES' then "Field_2" = 'True'
else "Field_2" = 'False'
end

but it returns NULL value.

I've tried with 't'and 'true' as well but it didn't work.

Is there anyone who knows how to solve it?


texts: Configuring custom network in the Directions dijit?

texts:
I'm making a CMV app and trying to make the directions widget work with my own network dataset. I replaced ESRI's street route network service with my own service, so viewer.js looks something like this:

directions: {
    include: true,
    id: 'directions',
    type: 'titlePane',
    path: 'gis/dijit/Directions',
    title: 'Directions',
    open: false,
    position: 7,
    options: {
        map: true,
        mapRightClickMenu: true,
        options: {
            routeTaskUrl: 'https://domain.com/arcgis/rest/services/my_network/NAServer/Route',
            routeParams: {
                directionsLanguage: 'en-US',
                directionsLengthUnits: units.MILES
            }
        }
    }

The app is not using my custom network. I get no results when using the widget--any suggestions as I troubleshoot?

I know the network service is working properly, as I got it to work in ESRI's stand-alone example here.


texts: PyQGIS. Field calculator based on the rank of another field

texts:
I am trying to make a script to calculate a field "rank" [RANK] based on the feature rank in another field "population" [POP]. My current script can found below. In this script, I created a new field "RANK".
The idea is that the city with the highest number of people will have the highest number in the field [RANK], namely 10. The city with the lowest number of inhabitants will have the lowest number, namely 1. So the idea is a kind of sorting number.

fn = path to the file

layer=iface.addVectorLayer(fn, '', 'ogr')
features=layer.getFeatures()
caps=layer.dataProvider().capabilities()

if caps&QgsVectorDataProvider.AddAttributes:
       res=layer.dataProvider().addAttributes([QgsField("RANK", QVariant.Double)])
   layer.updateFields()


texts: Any way to use the range slider in ArcGIS Pro with a floating point raster?

texts:
So I've created these classified floating point rasters showing specific heights. However, I need to enable range properties for them so that they can function with the range slider in ArcGIS Pro.

I am aware that this feature is only available for feature layers. I have tried converting the values to integers then using raster to polygon, however, that changed the values of many of them so the appearance was different (for example 3.9999 was now 3).

Is there a way then that I could directly transfer the rasters that I have into a different type of feature layer that could be range enabled? Maybe TIN or something else?

Their current appearance values must remain the way they are so I'm kind of in a bind.


texts: No difference between distances measured in projected and non-projected data

texts:
I'm performing measurements of distances in Km between a lot of localities by using the "Point Distance" in ArcGis 9.3. I realized that the results obtained after using projected coordinate system are the same of those obtained without projecting data. I change projected coordinate system in "Data frame properties". Where am I wrong?
I am a newbe in Gis and in ArcGis too. Any advices will be invaluable.
Thanks in advance.


texts: How to read *_LOW_HN / *_HIGH_HN fields in NY Street Centerline Data?

texts:
I've been looking at the NYC Street Centerline (CSCL), in particular number ranges for individual street segments.

Quite a few look like this:

           ST_NAME L_LOW_HN L_HIGH_HN R_LOW_HN R_HIGH_HN
0               28  215-001   215-027  215-000   215-026
1           QUEENS  120-011   120-011        0         0
2               58   51-001    51-099   51-000    51-098
3           FOREST   68-049    68-099   68-050    68-098
4               69  110-001   111-099  110-000   111-098
5              105  107-001   108-099  107-000   108-098
6           BEACON   14-025    14-099   14-024    14-098
7               49  108-001   110-099  108-000   110-098

an I have a trouble to wrap my head around these values.

How do I read for example 08-001 or 14-025? How do these correspond to odd and even numbers?

Can someone explain this on an example?


texts: Iterate zonal statistics table model builder - add field and name each row by raster

texts:
I have a model builder model to run zonal statistics as a table on multiple rasters and merge the output table into one table

However, I would like to add a field to the merged output table that identifies each each record by it's raster name.

Since I'm new to figuring this out, can someone tell me what steps I need to add in order to make that happen?

See image attached.

Thanks,


texts: Export file _Google Earth Engine_Null parameter passed to FileFilter get user data function

texts:
When I export Vector Shapefile or Raster file from GEE and try to open it on another software (ERDAS/QGIS) I am given the following error Null parameter passed to FileFilter get user data function.
Does anyone know what this means? I am conducting a flood analysis in GEE. The code for the export files is as follows:
// Export flooded area as TIFF file
Export.image.toDrive({
  image: flooded,
  description: 'Flood_extent_raster',
  fileNamePrefix: 'flooded',
  region: aoi,
  maxPixels: 1e10
});
// Export flooded area as shapefile
// Convert flood raster to polygons
var flooded_vec = flooded.reduceToVectors({
  scale: 10,
  geometryType:'polygon',
  geometry: aoi,
  eightConnected: false,
  bestEffort:true,
  tileScale:2,
});
// Export flood polygons as shape-file
Export.table.toDrive({
  collection:flooded_vec,
  description:'Flood_extent_vector',
  fileFormat:'SHP',
  fileNamePrefix:'flooded_vec'
});


texts: Processing a large number of raster files in QGIS

texts:
I have a large number (100,000+) of raster files I need to process in QGIS (exporting to tiff, clipping to polygon layers, etc), for a total of approx 500GB. Can anyone advise if there's some way of handling these within QGIS (a raster GDB?)? Can to stretch into code if needs be, but would prefer a simple solution using embedded tools / plugins if possible. I've used batch processing tools but some of the larger sub folders have 9000 files and I'm having problems loading them into the batch processing window.

Currently height (DSM) .asc files, UK coverage but patchy. Would prefer to work within QGIS as I'm more familiar with it, but would consider any open source / free toolset.


texts: Is there a way to export the style, along with a shapefile, to a kml in QGIS?

texts: I see there were plugins in the past that may have done this such as qgis2kml and qgis2google, but they no longer exist.


texts: Keep customized colors with geom_polygon()

texts:
I have plotted a map of France with customized colors with map():

map(database="france", fill=TRUE, col=color)

Results:

(don't mind the two errors)

and then I want to put this map with a projection of Lambert using ggplot2 package :

head(fortify(carte_france))
ggplot(map('france'), aes(long, lat, group = group)) +
  geom_polygon(col = 1) +
  coord_map(projection = "lambert",
            parameters = c(lat0 = 41.366005 , lat1 = 51.097523))

However it gives me this :

I know that the dark colours comes from the arg fill from geom_polygon(), but is there a way to tell the function geom_polygon() to not use the argument fill or to keep the colors I have put before?

Vector of colours:

color <- c("yellowgreen", "yellowgreen", "yellowgreen", "yellowgreen",
"gold1", "sienna3", "yellowgreen", "yellowgreen", "sienna3",
"gold1", "gold1", "gold1", "sienna3", "sienna3", "gold1", "dodgerblue2",
"dodgerblue2", "dodgerblue2", "gold1", "dodgerblue2", "sienna3",
"dodgerblue2", "burlywood2", "dodgerblue2", "palevioletred4",
"dodgerblue2", "palevioletred4", "dodgerblue2", "palevioletred4",
"palevioletred4", "gold1", "palevioletred4", "gold1", "darkorchid4",
"gold1", "darkorchid4", "palevioletred4", "green4", "cyan4",
"gold1", "palevioletred4", "burlywood2", "green4", "green4",
"cyan4", "darkorchid4", "green4", "burlywood2", "palevioletred4",
"burlywood2", "green4", "green4", "palevioletred4", "palevioletred4",
"palevioletred4", "green4", "burlywood2", "olivedrab2", "green4",
"olivedrab2", "darkorchid4", "darkorchid4", "cyan4", "darkorchid4",
"cyan4", "cyan4", "olivedrab2", "olivedrab2", "olivedrab2", "cyan4",
"cyan4", "olivedrab2", "cyan4", "olivedrab2", "olivedrab2", "cyan4",
"cyan4", "olivedrab2", "olivedrab2", "olivedrab2", "cyan4", "cyan4",
"cyan4", "cyan4", "orange2", "orangered1", "orangered1", "orangered1",
"olivedrab2", "orange2", "olivedrab2", "orangered1", "orange2",
"orangered1", "orange2", "orange2", "orangered1", "orangered1",
"orangered1", "orange2", "orangered1", "orange2", "orangered1",
"olivedrab2", "orangered1", "orangered1", "orangered1", "orangered1",
"orange2", "orange2", "orange2", "chartreuse3", "orangered1",
"chartreuse3")


texts: Opening Raster stored in Esri Geodatabase using QGIS

texts:
Is there a way to open a raster that is stored in ESRI V10.0 Geodatabase in QGIS V2.4?

I can open the vector layer, but there does not seems to be a way to open a raster.


texts: Is there a difference between the two sets of QGIS terrain analysis tools?

texts: Is there a difference between Raster/Terrain Anaylsis and the same features in Raster/Analysis/DEM (Terrain Analysis)? The two appear to do the same types of computations; are they the same or is there a difference and if so, what?


texts: Converting shapefile to raster: Character attributes

texts:
I've got a shapefile that I would like to convert into a raster, but the trick is that the attribute I want a raster version of is character data. I'm familiar with geoprocessing using GDAL (in the command line) and R, so solutions using either of those two are preferred. However, I have access to QGIS if necessary, and even ArcMap.

I'm using this shapefile. It contains polygons of municipalities in part of Brazil. The attribute NAME_2 is the name of the municipality, and I want a raster version of that.

I know that you can't make a raster with characters, so I created a numeric version of NAME_2 in R by doing

library("rgdal")
munisCrop.shape = readOGR(dsn=getwd(), layer="munisCrop")      # Read shapefile
munisCrop.shape$NAME_2_NUM = as.numeric(munisCrop.shape$NAME_2)  # Make new attribute
writeOGR(obj = munisCrop.shape, dsn = getwd(), layer = "munisCrop", driver = "ESRI Shapefile") # Save new version of shapefile

Values in munisCrop.shape$NAME_2_NUM range from 1 to 787.

I was hoping to then use gdal_rasterize on my command line (GDAL version 1.10.1 on Mac OS X 10.9.2) to make a raster version of the shapefile. Here's the code I was trying (separate lines for clarity):

gdal_rasterize -of GTiff -ot UInt32 -a NAME_2_NUM \
               -tr 0.0002777778 0.0002777778 -ot BYTE \
               -co COMPRESS=LZW -l munisCrop \
               munisCrop.shp munisCrop.tif

This completes easily enough. However, when viewing munisCrop.tif in QGIS, I see that the data range from 0 to 255. This is confirmed in R with

library("raster")
munisCrop.raster = raster("munisCrop.tif")
sort(unique(getValues(munisCrop.raster)))

0 makes sense, as NA values got converted to that. However, the maximum value of 255 is no good; remember, the maximum value of that attribute in the shapefile is 787.

I think I understand why this is happening but can somebody really explain it to me and suggest a workaround?

I've tried various options in gdal_rasterize instead of -ot UInt32 (for example, BYTE), and the result is always the same. I've also tried rasterizing it in R:

munisCrop.blankraster = raster(ext = extent(munisCrop.shp),
                               resolution = 0.0002777778,
                               crs = projection(munisCrop.shp)
                               )
munisCrop.Rraster = rasterize(x = munisCrop.shp,
                              y = munisCrop.blankraster,
                              field = "NAME_2_NUM"
                              )

But that seems like it takes forever. gdal_rasterize is way faster (a matter of seconds instead of... indefinite), so I'd love a solution using that.


texts:
I am trying to reproject some ascii files and convert them into shapefiles. I am new to python but looking at some code online and help files I was able to write the following code to run on Python console in QGIS.

import processing
for i in xrange(1, 191):
    filename1= "C:\\nsdbg\\GIS data\\watersheds\\ascii\\{}.asc".format(i)
    filename2= "C:\\nsdbg\\GIS data\\watersheds\\tifs\\basin_{}.tif".format(i)
    filename3= "C:\\nsdbg\\GIS data\\watersheds\\shapefiles\\basin_{}.shp".format(i)
    #reproject raster
    processing.runalg("gdalogr:warpreproject",filename1,
     "+proj=laea +lat_0=55 +lon_0=20 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs",
     "EPSG:3035",0,0,"",5,filename2)
    #raster into shapefile
    processing.runalg("gdalogr:polygonize", filename1,"DN",filename3)

As far as I can tell the reprojection is OK but when creating the shapefile I get 2 per ascii in different parts of the map. The one on the left (in the image below) with the right shape but on the wrong part of the map, and the one on the right in the right place but with a square around it (they are actually quite far apart in the map, I just collated the images).

If I use "polygonize" in the GUI instead of the python console this duplication does not happen and the location is correct. If I use it directly with the ascii file I get the right shape, however if I use it with the tiff file created by this code I do get the square around it (but not the duplication and it is in the right location).


texts: How to prevent QGIS automatically moves panel positions after the program restarts

texts:
I want to keep the user interface consistent. But each time I launch QGIS, the panel positions moves a little bit.
Below are three screenshots showing the first, second, and third time I launch QGIS. The panel positions moves downwards a little bit each time.
How can I fix it? I'm using QGIS 3.10 on windows 10 Home eidtion version 1809.


texts: Can't execute a simple psycopg2 script twice when running from ArcGIS Desktop as a script tool

texts:
I'm building a geoprocessing service that queries an SDE database (running on postgres) and returns results. I do a lot of development by calling the script from a toolbox in ArcGIS desktop 10.3.

For some reason when I call the script below in ArcGIS desktop, it executes successfully the first time, then fails every time after that with 'NoneType' object is not callable
on line 12 as I try to iterate over inAllCornerPointsList. If I close ArcMap and call it, it works the first time, then fails.

import psycopg2

conn = psycopg2.connect(host='localhost', database ='sde', user='sde', password='somepassword')
curs = conn.cursor()

sqlStatement = """SELECT * FROM sde.gly_ula_corners"""

curs.execute(sqlStatement)
inAllCornerPointsList = curs.fetchall()

for cornerPoint in inAllCornerPointsList:
    print cornerPoint

curs.close()
conn.close()

print 'Complete'

To further confuse things, if I select * from a different table (sde.countyboundaries), my script executes successfully every time.

The geometry type for sde.gly_ula_corners is shape st_point. The geometry type for sde.countyboundaries is shape st_geometry. This is the only difference that jumps out at me right away-- one is a point fc, and one is a polygon.

Any idea what might be causing this? Is there some weird memory environment at work? Calling the script multiple times from the command line or pythonwin does not fail.


texts: Use OSGeo4W shell command from python 2.7 script

texts:
I would like to convert jp2 files (sentinel-2) to geotiff format using python 2.7 on win 10. I saw some solutions using python and GDAL but they did not work for me. The only thing that did work was to open OSGeo4W shell, go to the desired directory using cd directory (from the shell) and then run (in the shell) for %i in (*jp2) do gdal_translate -of GTiff -co TFW=YES %i %~ni_tiff.tif. The problem is that I need to perform it on multiple folders and also there are some processes I need to perform on each file before and after the conversion, the conversion from jp2 to geotiff is just a small part of the process. So I would like to know if it is possible to call OSGeo4W shell command from within a python script. My desired results would be a python script that at some part of the script runs a command on the OSGeo shell, each time for a different folder (I can store the folders in a list if necessary) and continues to process the products of the OSGeo4W command shell.

I read solutions regarding command line (cmd) using os.system or subprocees but none of them referred to executing commands in a different platform than win command line.

The script according to my logic would be something like that:

import ...

#src is a list of directories that holds directories with jp2 files
src = ["c:\jp2_files_1", "c:\jp2_files_2","c:\jp2_files_3","c:\jp2_files_4"]

for i in (src):
    os.chdir(i)
    SomethingI'mMIssing('for %i in (*jp2) do gdal_translate -of GTiff -co TFW=YES
    %i %~ni_tiff.tif')

continue with python processing


texts: Meaning of NULL when using condition in gdal_calc

texts:
I am confused with the results of using a condition in gdal_calc.
Initially I thought that gdal_calc was setting pixels to nodata where a condition is not met.
In the documentation (http://www.gdal.org/gdal_calc.html) it is said that the following command sets values of zero and below to null.

gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>0)" --NoDataValue=0

If I try this command:

gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>0)" --NoDataValue=-9999

The pixels where the condition is not met are set to zero.

Should I understand that the term null used in the documentation actually refers to pixels with a value of zero?


texts: Convert exported ArcGIS .dwg to 3D rendering on AutoCAD 2013?

texts:
We're working on a bathymetry mapping project of our reservoir, digitized as it was built in 1977. I have a polyline .shp for contours and polygon .shp for area between contours with fields including depth, elevation, area for polygons and length for polylines. Unfortunately, we don't have Spatial or 3D Analyst. For visual presentation and ideally calculate max capacity, I believe the option we're left with is exporting as a .dwg to our Engineer's AutoCAD 2013.

Questions:

Is this our best option?
If so, what method should be used to convert to 3D rendering and what attribute data must be present to make it work?
Also, can z-values be stored in the 'elevation' field?


texts: Why is my Polygon (WGS84) in SpatiaLite 4 not displayed by QGIS?

texts:
I made a polygon in a SpatiaLite 4 database as follows:

INSERT INTO test_pg (id, name, geom) VALUES (3, 'test', ST_PolyFromText('POLYGON((8.556824 47.185269, 8.611009 47.185269, 8.611009 47.218568, 8.556824 47.218568, 8.556824 47.185269))', 4326));

I assume that the order of coordinates is always LON LAT - is that correct?

When I do SELECT ST_AsText(geom) FROM test_pg it correctly outputs the polygon coordinates.

Now in QGIS, I add the SpatiaLite layer but the polygon does not get displayed (even after clicking on this zoom to layer thing). Geometry seems to be valid, I checked it with the QGIS tool. I also found out that when moving the cursor on the map in the east direction, 47.xxx (LAT in this case) gets increased, when moving the cursor in the north direction, 8.xxx gets increased. This should be reversed, shouldn't it?

Where are my errors?

EDIT: I also INSERTED a polygon with LON LAT coordinates reversed. This polygon actually gets displayed in QGIS, but again, 47.xxx is displayed as LON and 8.xxx as LAT where it should be the opposite actually.

EDIT2: Consider this question as closed. For some reasons, it suddenly works. Don't ask me why, maybe restarting QGIS several times helped...


texts: I want to merge two polygon shapefiles into one without losing the attributes using ArcMap 10.1. To be exact I am creating a geological map and have some geologic formations in one shapefile and some other formations in another shapefile.


texts:
We are seeking an open source solution for dynamic (on the fly) IDW or Kriging points interpolation.

Basically we are looking for something similar to ArcGIS Spatial Analysis Service.

Interpolation can be performed in the browser or on the server side.


texts:
Attempting to gather MODIS aerosol data for a number of different AERONET stations in the UK to plot against one another. However, due to the number of files (upwards of 2000), QGIS and ArcMap crash.

I have attempted to find a solution within RSGISLib, whereby I extract a raster to .csv via a defined vector but keep getting the following:

  ERROR 1: Attempt to read shape with feature id (1) out of available
  range. Segmentation fault

I have seen a method to extract to a .csv using r, however, I would prefer for this to be done within Python. Is there a way?

import rsgislib
from rsgislib import zonalstats
import glob

inputimage = '/home/student/Desktop/Ben_Folder/Test2/test.kea'
inputvector = '/home/student/Desktop/Ben_Folder/AERONET_POINTS/Chilbolton.shp'
outputtxt = 'test.csv'
useBandNames = False
zonalstats.pointValue2TXT(inputimage, inputvector, outputtxt, useBandNames)


texts: How to get data from TileMill tileLayer properties in Mapbox for click/hover

texts:
I've successfully used Leaflet to load a geoJSON polygon-shapefile with data, and added click and hover interactivity to the layer. For a large dataset, I switched to mapbox.js 1.6.1, used TileMill to export a tiles, used infostream's mbtiles-php (https://github.com/infostreams/mbtiles-php) to serve the data.

I'm very confused, now, how to use mapbox.js to access show the polygon attributes in the tile layer? I have successfully loaded the tiles onto a page using mapbox.js - but don't know where to start to pull in the attribute table.

I did not find an appropriate Mapbox example on the example page, and the Mapbox support links are all directed to their contact page. Did I misunderstand how the purpose of Mapbox/TileMill?


texts: Arc Hydro toolbox fails to import to ArcGIS Pro

texts:
Has anyone had success importing the Arc Hydro toolbox or importing the archydro module in a python script?
tbx = arcpy.ImportToolbox(r"C:\Program Files\ArcGIS\Pro\Resources\ArcToolBox\toolboxes\Arc_Hydro_Tools_Pro.tbx")

---------------------------------------------------------------------------
SyntaxError                               Traceback (most recent call last)
In  [25]:
Line 1:     tbx = arcpy.ImportToolbox(r"C:\Program Files\ArcGIS\Pro\Resources\ArcToolBox\toolboxes\Arc_Hydro_Tools_Pro.tbx")

File C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\__init__.py, in ImportToolbox:
Line 165:   return import_toolbox(input_file, module_name)

File C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\toolbox_code.py, in import_toolbox:
Line 491:   mymodule = generate_toolbox_module(toolbox, None, False, False, False, module_name, use_alt_alias)

File C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\toolbox_code.py, in generate_toolbox_module:
Line 466:   'exec')

SyntaxError: EOL while scanning string literal (Arc_Hydro_Tools_Pro.tbx, line 6230)
---------------------------------------------------------------------------


texts: How to speed up processing large dataset

texts:
I have developed a processing tool in QGIS to quality check the data in 26 layers. The tool goes to each layer and checks so many conditions. when I run the tool for small data (less than 500 features in each layer) the tool is quite fast, only takes 3-4 minutes to run. But when I test this to large data (around 10000 Features on each layer) the tool takes forever to run (last test run 27 hours)
The tool currently accepts both GeoJSON and TAB files as input layers and the layers need to be opened in QGIS.
Can someone guide me on how can I increase the speed of this tool?
I cannot attach the full code as it is 9000 lines but the following show the structure of the code. The code checks every feature in every layer against some business rules
def initAlgorithm(self, config):

    parameter=QgsProcessingParameterMultipleLayers(
            'qgs_layers',
            self.tr('Input layers')
            )

    parameter.setMinimumNumberInputs(25)
    self.addParameter(parameter)

    parameter.setFlags(parameter.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
    self.addParameter(parameter)

    self.addParameter(
        QgsProcessingParameterFolderDestination(
            'Output_folder',
            self.tr('Output Folder'),

        ))

def processAlgorithm(self, parameters, context, feedback):

    for layer in parameters['qgs_layers']:
        layer_name=layer.name()

        try:
           exec('self.' + str(layer_name).replace(' ','_')+'_layer_obj=layer')
        except:
            pass

    self.QA_layer_1(feedback)

    self.QA_layer_2(feedback)

def QA_layer_1(self,feedback):
    features=self.structurepoint_layer_obj.getFeatures(QgsFeatureRequest().setFilterExpression (filter))
    for f in features:
        if feedback.isCanceled():
            sys.exit()

        if len(f['id'])!=18 and f['structure_type']!='Pole':
            self.error_list.append([table_name,f['id'],'Invalid ID',f['id'],'The StructurePoint ID should be 18 Character','point',f.geometry()])
            new_error=CriticalError()
            new_error.update({
             'table_name':'structurepoint','Feature ID':f['id'],'Error Type':'Invalid ID','value':f['id'],
             'description':'The StructurePoint ID should be 18 Character','geometry_type':'point','geometry':f.geometry()
            })

            self.error_dict_list.append(new_error)
            feedback.pushInfo('\t{} current value is ({})'.format(str(new_error['description']),str(new_error['value'])))

def QA_layer_2(self,feedback):
    pass


texts: Using a Select By Location in Python and Calculating a field based on what State Plane a point falls in

texts:
Here's what I have so far.  Do I need to make another feature layer based on the selection?

import arcpy

arcpy.MakeFeatureLayer_management(r"PATH TO THE POINT", "point_lyr")
arcpy.MakeFeatureLayer_management(r"PATH OF THE STATEPLANE", "stPlane_lyr")

arcpy.SelectLayerByLocation_management("point_lyr", "INTERSECT", "stPlane_lyr")

Here's the updated Script for each state plane zone

 #Feature layers
 Point_Layer = "Point_Layer"
 CP_Layer = "CP_Layer"
 IntersectLayerCON = "IntersectLayerCON"

 # Process: Make Feature Layer for the SMON Points
 arcpy.MakeFeatureLayer_management(SMONPoint, Point_Layer, "", "", "FID FID VISIBLE     NONE;Shape Shape VISIBLE NONE;link_fld link_fld VISIBLE NONE;State State VISIBLE NONE")
 # Process: Make Feature Layer for the Control Points
 arcpy.MakeFeatureLayer_management(ControlPoint, CP_Layer, "", "", "FID FID VISIBLE NONE;Shape Shape VISIBLE NONE;link_fld link_fld VISIBLE NONE;State State VISIBLE NONE")
 # Process: Make Feature Layer for Intersection Shapefiles
 arcpy.MakeFeatureLayer_management(IntersectLayer, IntersectLayerCON, "", "", "FID FID     VISIBLE NONE;Shape Shape VISIBLE NONE;OBJECTID OBJECTID VISIBLE NONE;ZONE ZONE VISIBLE NONE;ZONENAME ZONENAME VISIBLE NONE;FIPSZONE FIPSZONE VISIBLE NONE;SQMI SQMI VISIBLE NONE;COLORMAP COLORMAP VISIBLE NONE;Shape_Leng Shape_Leng VISIBLE NONE;Shape_Area Shape_Area VISIBLE NONE")

 #This will calculate the fields labeled the correct STATE PLANE zone
 # Process: Select Layer By Location for SMON
 arcpy.SelectLayerByLocation_management(Point_Layer, "INTERSECT", IntersectLayerCON, "", "NEW_SELECTION")
 # Process: Calculate Field for SMON
 arcpy.CalculateField_management(Point_Layer, "StatePlane", "\"Colorado North\"", "VB", "")
 arcpy.CalculateField_management(Point_Layer, "State", "\"Colorado\"", "VB", "")

It will just calculate the fields for state and state plane.


texts: Scale and density values in density map appear wrong when in projected coordinate system?

texts:
I'm making a series of maps illustrating different aspects (road density, location of dams, elevation etc...) of a single watershed in California. In most of my maps, I've had the projection set to the geographic coordinate system, NAD_1983 (I know it's better to use a projected coordinate system but some of my data layers for some of my maps just don't line up properly in the projected coordinate system and I don't have enough time or experience with GIS to figure out that problem).

However, I'm trying to make a map of road density, using the line density tool, and because I want my density units to be in km per square kilometer, I need to put my 'roads' layer in a projected coordinate system (otherwise, I can only use square map units). If I change all of my layers to a projected coordinate system (NAD_1983_California_Teale_Albers), the map itself looks great BUT the scale bar is totally off (in the geographic coordinate system, the scale is about 2in=20miles versus 2in=0.00007miles for the projected coordinate system). Everything I can find online talks about the scale bar being off when using a Geographic Coordinate system rather than a projection- which I understand-, but not when using a projected coordinate system instead of a geographic one which is my problem! I think the issue with the scale is also impacting my density maps because my density values are completely off (they should be maxing out at about 10km of road per square km and instead I'm getting values as high as 80km/km2).

Does anyone have advice either on how to fix my scale bar and length calculations when in the projection or even how to calculate density in metric units rather than map units in a geographic coordinate system?

I am using one of the most up to date version of ArcMap.


texts: Polygon from QgsPoints - QGIS 3.8 python

texts:
I'm trying to create a polygon layer from 4 points that I have already defined. There issue here is that I'm always getting the following error:

AttributeError: type object 'QgsGeometry' has no attribute 'fromPolygon'

These are the points (printed in console):

print (points)
<QgsPoint: Point (546016.50083752069622278 4760165.85990538075566292)>
<QgsPoint: Point (535433.25716883409768343 4749582.61623669415712357)>
<QgsPoint: Point (559882.5150161839555949 4746299.84572671726346016)>
<QgsPoint: Point (549299.27134749735705554 4735716.60205803066492081)>

And here is the code that produces the error above:

layer = QgsVectorLayer(str('Polygon?crs='+crs), 'polygon' , 'memory')
prov = layer.dataProvider()
points = [H1V1p,H1V2p,H2V1p,H2V2p]
feat = QgsFeature()
feat.setGeometry(QgsGeometry.fromPolygon([points]))
prov.addFeatures([feat])
layer.updateExtents()
QgsMapLayerRegistry.instance().addMapLayers([layer])


texts: exporting attribute tables of multiple rasters to excel

texts: I had a raster and I clipped it using a shapefile of polygons. Now I have about 5500 individual raster images. I have to export the attribute table of each raster to excel so that I can then append them and make a new dataset in excel. Can someone tell me how to do this?


texts: How to upload custom marker symbols to ArcGIS online?

texts:
There's a similar but not-quite-equivalent discussion here: Exporting symbology from ArcGIS Desktop to ArcGIS Online?

The people in that discussion, from 2011 to 2014, were using AGOL and trying to upload layers with symbology included.  I'm using a free individual account and trying to upload a custom image to symbolize a layer already displayed.  The interface seems to offer that option, but it doesn't respond as it should:

I've clicked on the "use an image" link indicated in the screencap.  In most interfaces this opens Finder, Windows Explorer, or similar.  In this interface all it does is open an input box.  I've tried pasting in the entire file path to Marker.png.  However, clicking on the blue circle containing the white plus sign, I consistently get:

  "Unable to load image."

The image is a simple .png with transparent background, 120 x 120 pixels.  I've tried Chrome and Firefox.  Anyone have any ideas?


texts: How to install PostGIS extension for PostgreSQL database under Ubuntu?

texts:
I am quite a newbie in using Postgres under Ubuntu.
I have created a PostgreSQL database and now want to install the PostGIS extension for it so that I can import shapefiles and use the PostGIS commands.
PostGIS is installed. I use the bash shell. Can somebody give me advice?


texts: How add a custom control panel in OpenLayers 3

texts: Is there an example that shows how to add a custom control panel in OpenLayers 3?


texts: String output as an input dataset variable

texts:
I have a model which creates topology with rules and validate itself.
But I got used to working with layers in TOC and have a script in Calculate value tool that with a help of Parse path tool make a string output with feature dataset path in.
Here is a script:

def GetPath(lyr):
    import os
    mxd=arcpy.mapping.MapDocument("CURRENT")
    target=arcpy.mapping.ListLayers(mxd,lyr)[0]
    return target.dataSource

And here is a part of model I want to be working.

Unfortunately, Create topology tool requires path not to some output value or variable, but probably direct file path. That is a reason for crashing of model due to the non-existance of dataset parameter.

Anyway is it possible to make an Input Feature dataset as a variable which records a path from string data? In my case it should be like %FDS%.


texts: How to create a map of demographics using Canada Census Data?

texts:
I understand this is an extremely vague question. I will start by telling you what I have:

I have Canada census boundary shapefiles and want to map data by dissemination area
I have Canada census .csv of populations and languages all organized by dissemination area
I am trying to join the .csv and the dissemination area shapefile table so I can show my boss what's up
I would like it to be organized kind of like what is shown here: https://www.realtor.ca/Residential/Single-Family/17597418/124-GADWALL-AV-BARRIE-Ontario-L4N8X3#v=de

Every time I try to make a query table it gives me this error - 000152 : Invalid input data. I have added new fields, tried changing the .csv, and now I am out of Ideas.

I got it from the StatCan website and I'm using ArcMap


texts:
Im trying to write some code (in C#, but any language will do) where the user defines and area, based on Min/Max Long/Lat's, and specifies a distance between points in meters, and the code returns a table of long/lats (WGS84/ESPG:4326) to create a grid.

I wrote some code a while ago that I was under the impression did the job, but after a quick bit of analysis, doesnt quite measure up to the job.  Please see the code below:

decimal minX = -74.002747535706M;
decimal minY = 40.722282672831M;
decimal maxX = -73.98386478424M;
decimal maxY = 40.733470232685M;
decimal distance = 100;
decimal x = 0.0M;
decimal y = 0.0M;
decimal length = 0.0M;
decimal lengthIncrease = 0.0M;

decimal increasePercentage = 0.0M;
int i = 0;

decimal yIncrement = distance / (decimal)111111.111111;

DataTable dt = new DataTable();
dt.Columns.Add("Longitude" , Type.GetType("System.Decimal"));
dt.Columns.Add("Latitude" , Type.GetType("System.Decimal"));

y = minY;

while (y < maxY)
{
    while (x < maxX)
    {
        length = (decimal)111.325 * (decimal)Math.Cos((double)y * (double)0.0174532925199433);
        increasePercentage = (decimal)0.001 / length * 100;
        lengthIncrease = increasePercentage / 100 * distance;
        x += lengthIncrease;
        i++;
        DataRow dr = dt.NewRow();
        dr["Longitude"] = x;
        dr["Latitude"] = y;

        dt.Rows.Add(dr);
    }
    y += yIncrement;
    x = minX;
}

return dt;

The above works fairly well.  It returns long/lats seperated by 100m lat, but by 90m long.  When i change the distance to 25m, it seperates latitude by 25m, but longitude by 22.3m.

EDIT:  C Sharp Version Solution Based on yosukesabai Answer:

int radius = 6378137;
double pi = Math.PI;
double deg2Rad = pi / 180;
double dst2Lat = 360 / (2 * pi * radius);

DataTable dt = new DataTable();
dt.Columns.Add("ID" , typeof(int));
dt.Columns.Add("Longitude", Type.GetType("System.Decimal"));
dt.Columns.Add("Latitude", Type.GetType("System.Decimal"));

x = minX;
y = minY;
while (y < maxY)
{
    double rLat = radius * Math.Cos(deg2Rad * y);
    double dst2Lon = 360 / (2 * pi * rLat);
    while (x < maxX)
    {
        i++;
        x += dst2Lon * distance;
        DataRow dr = dt.NewRow();
        dr["ID"] = i;
        dr["Longitude"] = x;
        dr["Latitude"] = y;

        dt.Rows.Add(dr);
    }
    x = minX;
    y += dst2Lat * distance;

Any help here would be fantastically helpful!

Regards

AM


texts:
I am trying to generate a Cost Distance table using Spatialite.
I have several points in the table random_points for which I would like to compute the Delaunay triangulation and use the segments it produces as arcs.
Something like:

SELECT ST_DelaunayTriangulation(ST_Collect(Geometry)) FROM random_points;

And the reuse its output. How can this be achieved?

Alternatively, I have thought about generating the lines connecting each pair of random points (there will be millions of these, the total number of random_points is around 4000) and then to keep those that appear in the Delaunay triangulation. Something of that sort:

WITH dela AS (SELECT ST_DelaunayTriangulation(ST_Collect(Geometry), 1, 0) FROM random_points),
  g1 AS (SELECT mergename, geometry FROM random_points),
  g2 AS (SELECT mergename, geometry FROM random_points)
  SELECT g1.mergename, g2.mergename, AsText(MakeLine(g1.geometry, g2.geometry))  FROM g1, g2 WHERE MakeLine(g1.geometry, g2.geometry) IN dela;

Any idea how to tweak it if the syntax is right?


texts: Distortion at edges of planetscope images

texts: Does anyone know why there is distortion on the edges of some planetscope images?


texts: How to read coordinates from shapefiles and overlay on Google Maps?

texts:
I have shapefiles i.e. *.shp, *.dbf and *.shx files.

I need to read  geo coordinates from these files and overlay the area on Google Maps at the right location.

I believe that *.shp and *.dbf files contain coordinates.

I have also used gotools api to read ArcGIS files and am able to generate map.

Is it possible to do this?


texts:
I'm currently working with Landsat 8 OLI/TIRS Level-2 Data Products - Surface Reflectance and I would like to compare different landsat images from different places on earth. I was wondering if pixels in this images are comparable or I should re-scale or pre-process the image using the metadata of the file or something, or I can directly compare pixels of images across time and space.

According to the LaSRC Product guide, the unit is reflectance and the range goes from -2000 to 16000, but shouldn’t it takes values from 0 to 1 or something like that?

Would it be safe to re-scale this as it goes?

Or should I use the metadata file to check for maximums and minimums?

I’m sorry. I just sometimes find literature a bit confusing or inconsistent with the terminology, so I’m not sure.


texts:
I have a plugin that loads a Raster image and selects the correct base maps to use.
The Raster image load with no problem, but my base maps will not load.

I have create a qlr, which has 5 references. I can drop it into the Layer pane, so I know that the file is good.

What is wrong with my script?

# get a handle to the layer groups form the default project
rootLayers = QgsProject.instance().layerTreeRoot()

 # load an group the base layers
 vlayer = QgsVectorLayer(self.Settings.staticLayer,"static","ogr")

 # load and group the raster
 rlayer = QgsRasterLayer(self.rasterFile.fileName(),"analysis")

 QgsMapLayerRegistry.instance().addMapLayer(vlayer)
 rootLayers.findGroup("Static").addLayer(vlayer)
 QgsMapLayerRegistry.instance().addMapLayer(rlayer)
 rootLayers.findGroup("Analysis").addLayer(rlayer)


texts:
I have a GeoJSON file called mygeojson.json and I want to add it as a layer in OpenLayers 3 on top of an openstreetmap layer. So far I could display the openstreetmap world including zoom etc. but for some reason I can't get the mygeojson.json on it.

The geojson contains many polygons and looks like this:

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

"features": [
      { "type": "Feature", "properties": { "DN": 2 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 13.559093915055664, 52.545214330050563 ], [ 13.559633429050496, 52.545205649772548 ], [ 13.559633415380715, 52.545214636296755 ], [ 13.559093915055664, 52.545214330050563 ] ] ] } }
]
}

my main.html:

<!doctype html>

    <link rel='stylesheet' href='http://ol3js.org/en/master/css/ol.css'>

    <title>OpenLayers 3 example</title>

My Map

I also tried removing the projection info but no use.


texts:
I have installed Postgresql 9.3.5 and Postgis 2.1.2 on a 16Gb Ubuntu 14.04 system.

(1) If I load all the states TIGER data, using the scripts provided in Postgis, how big will the database be for the entire US?

I have currently installed just the TIGER data for California (CA).

(2) It takes about a minute to geocode an address.  I have read the performance hints suggested at How to improve poor performance times for PostGIS 2.1 Tiger Geocoder?, and have tried setting the shared memory buffer to a much larger value.  I haven't been able to improve the query performance time.

Do you have suggestions?


texts:
Will there be available stable version of QGIS 1.7.4 in Ubuntu 12.04?

Version 1.9 is working fine at the moment, but as it is unstable currently,  don't want to risk it with my work.


texts:
My initial code read an image, apply a median filter to it, and save the smoothed image. It looked like this :

#Image reading
rasterpath = r"C:/Users/SIEXT5646/Desktop/Donnes/2018_gris.tif"
ds = gdal.Open(rasterpath)
band1 = np.array(ds.GetRasterBand(1).ReadAsArray())

#Median filter computation :
im_med = ndimage.median_filter(band1,size=20)

#Export of smoothed image

However I am working with very large image (orthophotos of infrastructure, such as tunnels), so I would like to add a process bar to my code. I used the idea developed in the link here : Displaying progress bar when using QGIS processing tool. So now my code looks like this now :

# Progression bar
qgis.utils.iface.messageBar().clearWidgets()
progressMessageBar = qgis.utils.iface.messageBar()
progressbar = QProgressBar()
progressMessageBar.pushWidget(progressbar)

# Processing feedback
def progress_changed(progress):
#print(progress)
progressbar.setValue(progress)

f = QgsProcessingFeedback()
f.progressChanged.connect(progress_changed)

#Read image

#Median filter computation, with the addition require for the feedback
processing.runAndLoadResults(im_med = ndimage.median_filter(band1,size=20), feedback=f)

#Export smoothed image

The error I'm getting is TypeError: runAndLoadResults() got an unexpected keyword argument 'im_med', so I'm guessing the problem is coming from the fact that I define a new variable Inside the processing.runAndLoadResults.

Does anyone know how to solve this ?


texts: Syntax error running script in PythonWin after success in ArcPy

texts:
I am building a script to automate a water outage map(I'm new to Python). In the ArcGIS Desktop 10.1 window the script runs fine. However in Pythonwin and IDLE I get syntax errors(In Pythonwin reads: Failed to run script-syntax error-Invalid syntax). The goal is to use task scheduler to run the script outside of ArcGIS.

import arcpy
mxd = arcpy.mapping.MapDocument(r"X:\Mikes_Workspaces\Outage\Outage.mxd")
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
# Script arguments
OutCurrent = arcpy.GetParameterAsText(0)
if OutCurrent == '#' or not OutCurrent:
    OutCurrent = "X:\\Geodatabases\\WebData\\Water_Service.gdb\\OutCurrent" # provide a default value if unspecified

# Local variables:
Service_Group = "Service_Group"
Update_ = "X:\\Mikes_Workspaces\\Outage\\Outage_Board.xls\\Update$"
Group_Out = "Service_Group"

# Process: Add Join
arcpy.AddJoin_management(Service_Group, "Group_", Update_, "Group_Out", "KEEP_COMMON")

# Process: Copy Features
arcpy.CopyFeatures_management(Group_Out, OutCurrent, "", "0", "0", "0")
    # Process: Symbology
arcpy.ApplySymbologyFromLayer_management("OutCurrent", "X:\Mikes_Workspaces\Online Shapefiles\Outage_today.lyr")
# Process: Remove Join
arcpy.RemoveJoin_management(Service_Group, "")
mxd.save()
arcpy.RefreshActiveView()


texts: Converting the distance in meters using object.distance function from shapely lib

texts:
This code allows to obtain the distance from a data point to each element of my geojson (it uses real coordinates, the ones used by google maps), and then I can know which is the closest point.
In the image, red point is my input coordinates, and the shape with border red is the nearest point according my logic.

In this location I am getting this value 0.0002648... in properties in my array of json features[x].properties.distance. I need convert that value in meters.
I want to make a condition according to the distance meters from the given point.
How can I do the conversion to pass this value to meters?
import json

from shapely.geometry import Point, mapping, shape
from shapely.ops import nearest_points

with open("g.geojson") as f:
    geojson = json.load(f)

ref_point = Point(-95.29756, 29.85736)

features = []
min_distance = None
min_index = None

for i, feature in enumerate(geojson["features"]):

    # Could be improved with Haversine distance, or distance in UTM coordinates for example
    p1, p2 = nearest_points(shape(feature["geometry"]), ref_point)
    distance = p1.distance(p2)

    if min_distance is None or distance < min_distance:
        min_distance = distance
        min_index = i

    feature["properties"] = {"index": i, "distance": distance}
    features.append(feature)
    features.append(
        {
            "type": "Feature",
            "properties": {"index": f"closest_point_{i}", "distance": distance},
            "geometry": mapping(p1),
        }
    )

features = sorted(features, key=lambda x: x["properties"]["distance"])

features.append(
    {
        "type": "Feature",
        "properties": {"index": "ref_point", "distance": 0},
        "geometry": mapping(ref_point),
    }
)


texts: Unable to add shapefiles via drag and drop - QGIS 2.4 OSx

texts:
I have a clean install of OSx Yosemite Beta 5, with homebrew install of QGIS 2.4 and all its basic dependencies.

-QGIS version    2.4.0-Chugiak
   -QGIS code revision  exported
   -Compiled against Qt 4.8.6
   -Running against Qt  4.8.6
   -Compiled against GDAL/OGR   1.11.1
   -Running against GDAL/OGR    1.11.1
   -Compiled against   GEOS 3.4.2-CAPI-1.8.2
   -Running against GEOS    3.4.2-CAPI-1.8.2 r3921
   -PostgreSQL Client Version   9.3.5
   -SpatiaLite Version  4.2.0
   -QWT   Version   6.1.1
   -PROJ.4 Version  480
   -QScintilla2 Version 2.8.3

Am able to load QGIS with no errors thrown, and can view shapefiles, connect to postgres, etc, but can only add new shapefiles if done through the 'Layer - Add Vector Layer' workflow.  Dragging and dropping into QGIS results in:

'Invalid Data Source:/.file/id=6571367..... is not a valid or
recognized data source

I am also unable to use SPIT to import shapefiles into postgres:

Problem inserting features from file:/xxxxx/xx.shp


texts: Missing zoom levels in MBTiles generated from GeoServer

texts:
I have added MBTiles plugin to my GeoServer in order to generate my maps in .mbtiles format, so I can read maps in offline mode.

My problem is that, when I preview my maps in GeoServer, it includes all zoom levels, but once I downloaded .mbtile file I found that it contains only zoom levels from 1 to 5. And when I changed the EPSG:4326 to EPSG:3857 it doesn't work but when I opened the .mbtiles file from dbBrowser for SQLite 3 I found that it contains another range of zoom levels: this time from 19 to 21.

The request is:

http://localhost:8080/geoserver/wms?request=GetMap&service=WMS&version=1.1.1&layers=mdintimaps:test&styles=&srs=EPSG%3A4326&bbox=-76.962,-76.735,120.066,69.01&width=768&height=568&format=mbtiles

Any help?


texts: How to create a custom data viewer on top of QGIS?

texts:
I was exploring the possibility of creating a simple Data Viewer based on QGIS.

I would like to create an simple application that will open QGIS up with a limited set of tools & menus. QGIS would open up with a given Project file and data (in shapefile format).

Preferably, I'll like to make a stand alone installation program, that will configure QGIS and copy the data like I have mentioned above.

Is something like this possible?


texts: Extract values from raster from list of text file of coordinates in r

texts:
I have several raster files that i am trying to extract values from that come from some lat longs found in a text file. I've been successful with raster::extract in the past but seem to be running into a new problem here. The extract function returns NAs despite the projections seeming to line up.

Here is an (almost) reproducible example that shows my issue.

#reproducable example for SO
#
library(leaflet)
library(rgeos)
library(rgdal)
library(raster)
testRast <- raster("FranceRasters/FRA_STU_EU_S_CLAY.rst")
# crs string for ETRS89
crs(testRast) <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

testRast
# class       : RasterLayer
# dimensions  : 1109, 1076, 1193284  (nrow, ncol, ncell)
# resolution  : 1000, 1000  (x, y)
# extent      : 3209000, 4285000, 2026000, 3135000  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
# data source : D:\GIS\STU_EU_Layers\FranceRasters\FRA_STU_EU_S_CLAY.rst
# names       : FRA_STU_EU_S_CLAY
# values      : 0, 61  (min, max)

#Small sample of lat ongs that need sampling
Coordinate.sample <-  structure(list(Longitude = c(7.62648877, 7.63023291, 7.37533534,
                                                   7.37484746, 7.59758433, 7.46470834, 7.46133708, 7.46632616, 7.45291197,
                                                   7.52007514), Latitude = c(48.56788166, 48.59294834, 48.65367653,
                                                                             48.65288703, 48.78729009, 48.71844903, 48.71878129, 48.71828592,
                                                                             48.71898041, 48.52164554)), .Names = c("Longitude", "Latitude"
                                                                             ), row.names = c(NA, 10L), class = "data.frame")

#convert to spatial points
Coordinate.sample <- SpatialPoints(Coordinate.sample)

projection(Coordinate.sample) <- projection(testRast)
#show projections line up
leaflet() %>% addTiles() %>% addRasterImage(testRast) %>% addCircleMarkers(data=Coordinate.sample)

#extraction leads to NAs
raster::extract(testRast,Coordinate.sample)
      [,1]
 [1,]   NA
 [2,]   NA
 [3,]   NA
 [4,]   NA
 [5,]   NA
 [6,]   NA
 [7,]   NA
 [8,]   NA
 [9,]   NA
[10,]   NA

As you can see the extraction has failed. Does anyone have any ideas why it has failed?

Here is an image of the leaflet() output to show that the coords and raster are in the right projection.


texts: Raster Calculator for monthly input rasters from postgis

texts:
I have a PostGIS table (soil_moist) with columns as

id (integer), rid (integer), fdate (date), rast(raster), layer(integer)

The table is updated every day for 3 layers (layer value as 1,2 and 3). Hence, there are 3 rows added every day for the same area. The table contains the record from 2009-01-01 up to 2016-12-31. My task is to calculate the monthly average raster for layer 1 from 2009-01.

So, I created table as soil_moist_avg and then tried to accomplish this using Python and QGIS. My code snippet is as follows:

# -*- coding: utf-8 -*-

from datetime import timedelta, date
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

config = {
    "dbname"  : "model",
    "user"    : "postgres",
    "host"    : "localhost",
    "password": "******"
}

start_date = date(2009, 1, 1)
end_date = date(2017, 1, 1)
delta = timedelta(1)

initial_date = start_date
layer_entry = []
num = 1
while start_date < end_date:
    connString_ = "PG: dbname=%(dbname)s host=%(host)s user=%(user)s password=%(password)s port=5432 mode=2 column=rast table=soil_moist where=%(where)s" % \
                      {"dbname": config["dbname"],
                       "host": config["host"],
                       "user": config["user"],
                       "password": config["password"],
                       "where": 'fdate=\'%s\' AND layer=1' % start_date
                       }
    layer = QgsRasterLayer(connString_, "%s_%s" % (str(start_date), str(num)))
    raster_entry = QgsRasterCalculatorEntry()
    raster_entry.ref = 'layer@%s' % num
    raster_entry.raster = layer
    raster_entry.bandNumber = 1
    layer_entry.append(raster_entry)
    num+=1
    if ((start_date + delta).month > start_date.month) or ((start_date + delta).year > initial_date.year):
        output_path = "D:/output/%s.tif" % start_date
        expression = "("
        for idx, layer_ in enumerate(layer_entry):
            expression += "layer@%s*" % str(idx+1)
        expression += "1)"
        expression += "/%s" % str(len(layer_entry))
        calc = QgsRasterCalculator(expression, output_path, "GTiff", layer.extent(), layer.width(), layer.height(), layer_entry)
        calc.processCalculation()
        # Reset Values
        num = 1
        layer_entry = []
    start_date+=delta

Here, I am trying to dump the raster as tif for the monthly average and later thought of dumping into PostGIS database. However, the layer returned from the connection string is empty (I checked with layer.width() and layer.height()).

What am I doing wrong here?

Is it because of the where clause?

Also, I would like to know if there is a better way to accomplish this task, to directly dump the output raster into PostGIS raster data type?


texts:
I've received a Arcgis layer containing streets, represented by polylines. When I identify a specific street, a large attribute table pops up:

OBJECTID
Shape
LINK_ID
ST_NAME
FEAT_ID
ST_LANGCD
NUM_STNMES
ST_NM_PREF
ST_TYP_BEF
ST_NM_BASE
ST_NM_SUFF
ST_TYP_AFT
ST_TYP_ATT
ADDR_TYPE
L_REFADDR
L_NREFADDR
L_ADDRSCH
L_ADDRFORM
R_REFADDR
R_NREFADDR
R_ADDRSCH
R_ADDRFORM
REF_IN_ID
NREF_IN_ID
N_SHAPEPNT
FUNC_CLASS
SPEED_CAT
FR_SPD_LIM
TO_SPD_LIM
TO_LANES
FROM_LANES
ENH_GEOM
LANE_CAT
DIVIDER
DIR_TRAVEL
L_AREA_ID
R_AREA_ID
L_POSTCODE
R_POSTCODE
L_NUMZONES
R_NUMZONES
NUM_AD_RNG
AR_AUTO
AR_BUS
AR_TAXIS
AR_CARPOOL
AR_PEDEST
AR_TRUCKS
AR_TRAFF
AR_DELIV
AR_EMERVEH
PAVED
PRIVATE
FRONTAGE
BRIDGE
TUNNEL
RAMP
TOLLWAY
POIACCESS
CONTRACC
ROUNDABOUT
INTERINTER
UNDEFTRAFF
FERRY_TYPE
MULTIDIGIT
MAXATTR
SPECTRFIG
INDESCRIB
MANOEUVRE
DIVIDERLEG
INPROCDATA
FULL_GEOM
URBAN
ROUTE_TYPE
DIRONSIGN
EXPLICATBL
NAMEONRDSN
POSTALNAME
STALENAME
VANITYNAME
JUNCTIONNM
EXITNAME
SCENIC_RT
SCENIC_NM
TO_X_LANES
FR_X_LANES
FOURWHLDR
COVERIND
PLOT_ROAD
REVERSIBLE
EXPR_LANE
CARPOOLRD
PHYS_LANES
VER_TRANS
PUB_ACCESS
LOW_MBLTY
PARK_AVAIL
PRIORITYRD
SPD_LM_SRC
TRANS_AREA
L_FIPS
R_FIPS
L_STATE
R_STATE
L_CITY
R_CITY
MEANDER
DRIVE_TIME
L_ROUTE
R_ROUTE
FROM_ELEV
TO_ELEV
OLD_LROUTE
OLD_RROUTE
BREADTH
WALK_TIME
WALK_PEN
NO_PARK_L
NO_PARK_R
TRAVCONS
LEFT_NOTE
RIGHT_NOTE
FT_PENALTY
TF_PENALTY
FUTURNRES
TUTURNRES
AR_MOTOR
EXPAND_INC
DCA
DISTANCE
Serv_Time
RSOwner
RSClass
RSDivision
RSRoadClass
RSPaveClass
Route_Num
Beltway
Service Left
Service Right
# of Passes Left
# of Passes Right
Left Priority
Right Priority
Left Range
Right Range
RS_Truck
RS_ServiceLeft
RS_ServiceRight
RS_TUTURNRES
RS_FUTURNRES
RS_MISC
RS_5T_TRAVCONS
RS_1T_TRAVCONS
RS_HOOD_L
RS_HOOD_R
Shape_Length

Although I suspect that there are some custom attributes in this table, the majority of the attributes in the table are listed here: http://webapub.sjrwmd.com/arcgis/rest/services/srvc/srRegMapExcnty/MapServer/17

Is there a way to find a description of these attributes? Searching for 'arcgis street attributes/fields/...' does not yield any solutions. Since I'm new to arcgis, I've no idea where to start searching. Searching for some specific attribute names did not yield any results. This external page http://sumo.dlr.de/wiki/Networks/Import/ArcView gave some insights, but only a few attributes are defined there.

Similarly, is there somewhere more information on a polyline? I need to create a graph, exactly like this: Generating graph edge list from polyline and point shapefiles?
Although I understand the proposed solution, it seems really error-prone to match points to end-points of edges based on their geographical coordinates. I don't understand why edges are not simply specified by (u,v), where u and v are identifiers of the starting and ending nodes. Two edges (i,j) and (u,v) are then connected if j=u. At least this is how it's done in Openstreet maps.


texts: How can I get the street names (in vector layer OSM, lines) that border parcels of land (in vector layer, polygons) on ArcMap?

texts: I have parcel data and street-level data (not in map), which I need to merge. My map contains the polygon-parcels and OSM streets - is there a way to capture the 4 or 5 streets that surround a parcel on ArcMap?


texts: Extract linestring from geometry collection

texts:
In QGIS, at SQL query intersection results in a geometrycollection column with points, linestrings and multilinestring. (first image)

I only need linestrings and multilinestring geometry to calculate the lenghts.

Using ST_CollectionìExtract I tried to select only linestring geometries but I obtain linestring empty for point geometry. (second image).


texts:
I want to create a new layer/shapefile, containing polylines that represent a connection to a grid. I have a layer of points and a layer of polylines (see picture below), based on this I want to connect each point to the perpendicularn (closest) polyline. These lines should then be stored in a separate shapefile.

Someone can help me out here? I am using ArcMap 10.4.1 - Basic Edition (Desktop version) (moving to standard edition in a couple of days)

As I work with a big dataset (3 mio points / 700 000 lines), I don't want to snap line per line.

EDIT: I will not be able to upgrade to the advanced edition and therefore not able to use the 'Near' tool


texts:
Is bulk conversion from csv to shapefile possible with ogr2ogr? I can do it file by file when I first create dbf file, then using vrt file I create shapefile. But using command

for /R %f in (*.csv) do ogr2ogr -t_srs EPSG:4326 -f "ESRI Shapefile" "%~dpnf.shp" "%f

I am able to get only the dbf and prj file, but not qpj, shp nor shx files. I have hundreds of csv files in one folder I would like to convert to shapefiles, but I'm on half way to achieve this.

Edit
For a single file I first created .dbf file: ogr2ogr -f "ESRI Shapefile" all8.dbf all8.csv

Then, I created all8.vrt file:

<OGRVRTDataSource>
  <OGRVRTLayer name="all8">
    <SrcDataSource>all8.csv</SrcDataSource>
    <SrcLayer>all8</SrcLayer>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>WGS84</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="Longitude" y="Latitude"/>
  </OGRVRTLayer>
</OGRVRTDataSource>

And finally run: ogr2ogr -f "ESRI Shapefile" all8 all8.vrt

I also tried to replace "all8" in .vrt file with "*" (to create multiple shapefiles using only one .vrt file), but that did'nt work.