texts:
I have a shapefile of all US States that I've imported into a table with a multipolygon geometry column. (SRID 4326)

I perform a spatial query based on lat/lon to return the state abbreviation. This normally works just fine. Today a query came in for Hawaii and the query returned nothing. So naturally I opened the table up in QGIS and queried the point's coordinates. Sure enough it returned the correct result. I clicked identify and it most certainly had the correct attributes.

The query I'm using is:

SELECT states.state_abbr FROM states WHERE st_contains(states.geom, ST_SetSRID(ST_MakePoint(-159.717,21.933), 4326))

Does anything look wrong with that query? I tried st_intersects as well.


texts: Fixing non-functional ArcMap which was OK day ago?

texts:
For about 1 month ArcMap 10.6 was working no issue, in conjunction with Intro to GIS online course. For the past 24 hrs.—approaching final project deadline—data files display, but with Nonsense Attribute Tables, No Data Manipulation function, etc.

I eventually noticed that despite saving all GIS data to file output:
C:\GIS\IntroductionToArcMap_Data\data,  ArcGIS was duplicating and filling with also unknown files at a default geodatabase path:
C:\Users\Owner\Documents\ArcGIS\Default.gdb

I tried changing path default (1: deleting all GIS files 2: setting Current and Scratch Workspace to outputs via Geoprocessing Environment Settings to paths instructed with initial download data instructions  3: Downloaded initial data (containing source, scratch & results .gdb files)into those paths.)
ArcMap worked for again for about 20min. Still sending out data to same default path—no save or export—just opening data in ArcMap from course download.

I am writing this before I attempt to delete & reinstall ArcMap & all GIS data, but with little optimism that it resolves the issue.

Can anyone offer actual insight into cause(s)/solution(s)?


texts: Creating TIN and contours in ArcGIS for Desktop?

texts:
I am upgrading from ArcView 3.1 to ArcMap 10.1. I have multiple LIDAR (bare earth) .txt files (x,y,z) that I would like to use to create a "TIN" (or is it Terrain Dataset in 10.1) and then contours?

I am unable to figure out how to do this in 10.1 and I have googled and read just about everything.

It was easy in ArcView 3.1 even though it had limitations (just a few steps). I have the full licensed 10.1 with all the 3D Analyst and variable extensions.

Can anyone point me in the right direction as I am currently working on a project?


texts:
I am currently working to convert our cities storm and sanitary sewers from 2D to 3D. Currently, we create profile views in CAD which means a lot of manual work anytime anything changes. We have a considerable number of manholes and mains mapped in GIS, and the ability to have dynamic data seems invaluable when it comes to profile work.

What I'm wondering is what the best way to go about this work is. I've been trying to use the point profile tool from the 3D analyst extension, but I'm unable to use it while showing our mains. I have downloaded the water utility network extensions, but I haven't been able to find any good documentation about how to actually use these or if I need to reformat all my existing data to use it.

I am lucky enough to have an enterprise license and should have access to any ESRI tool. We're running Arc 10.4. I know there's many different approaches so I'm very open to suggestions.


texts:
I am working on a web application with OpenLayers 6 currently, and it is my first foray into OpenLayers. I have managed to get the application to allow drawing and retrieval of images based on those bounding boxes. I understand the map has layers, and the Vector layer to an extent is what allows the drawing interactions / houses the features if I am not mistaken. What I do not get is the hierarchy of managing images through an Image layer and see no documentation.

Does each image have to be it's own layer (ImageLayer)?

I have my map initialized like this:

this.olmap = new Map({
  target: null,
  layers: [this.raster, this.vector],
  controls: [],
  view: new View({
    center: this.state.center,
    zoom: this.state.zoom
  })
})

and when features are drawn they are added via this.olmap.addInteraction(feature) which adds them to this.vector. Makes enough sense to get by.

When adding an image to the map via WMS I do so like:

const featureImage = new ImageLayer({
  source: new ImageWMS({
    ratio: 1,
    url: `https://myurl.com`,
    params: {
      LAYERS: id,
      TILED: true,
      serverType: 'geoserver',
      transition: 0,
    }
  })
})
featureImage.setOpacity(0.5)
this.olmap.addLayer(featureImage)

as that is the only way I could get to add an image successfully. But the resulting layers array continues to grow the more images are queried (obviously)

Now I did try

this.imageSource = new ImageWMS()
this.image = new ImageLayer({source : this.imageSource})

and add that to the initial layers of the map which works fine, but I cannot seem to add images to that subsequently via setting the source of this.image but seems like that would just override each iteration.

At the end of the day I suppose I could just add ImageLayer after ImageLayer to the map and manage them via id's in the array of layers, but that just seems dumb, and there should be a better solution, but I guess as you know starting out with OpenLayers is daunting. Looking for any advice on layer management that the API does not have, and how I could go about managing multiple ImageWMS() objects under one layer.


texts: Selecting points within a polygon in PostGIS

texts:
In PostGIS, we are trying to select the rows in a table of points (point_table) that are within a polygon, which is stored in a separate relation (poly_table). Having dug around for hours, we can't work out how to get our desired result, all the points within each polygon.

The closest we've gotten (i.e. no errors) is below, which returns 0 rows:

SELECT * FROM point_table as p
WHERE ST_Within(p.point::geometry,
(SELECT polygon FROM poly_table WHERE objectid=1)::geometry);

We did get a result of all FALSE boolean values from the query below, which we think is why the query above returns 0 rows...

SELECT ST_Within(p.point::geometry,
(SELECT polygon FROM poly_table WHERE objectid=1)::geometry)
FROM point_table as p;

We added these explicit type casts because any alternative threw an error, even though both the point and poly attributes are geography datatypes.


texts: Is there any plugin to Identity in QGIS?

texts:
Identity is one of ESRI ArcGIS Feature that clipping feature class by another feature class then union first feature class to the result of clip process.

The simple explanation is like this picture below.

Is there any plugin or feature like this in QGIS?
I'm using QGIS 2.14 Essen.


texts: Date formating in QGIS using expression?

texts:
My expression:

day(now()) || '/' || month(now()) || '/' || year(now())

returns current date as 15/4/2019, but I want to get the date in dd/mm/yyyy format i.e. 15/04/2019.

How can I get that using QGIS desktop version 2.18.12?


texts: How to generate a grid of horizontal equal-area rectangles in a Mercator projection?

texts:
I'm new to GIS so pardon my ignorance and bear with me. I'm generating a grid relative to a LonLat point. To make the cells all equal area, I'm creating it in 2163 projection. The problem is that when I reproject to 900913 and stick it on a Google Maps baselayer in OpenLayers, the whole grid is rotated clockwise. I understand this is a purely cosmetic issue, but it's a project requirement. Is it possible to change the map's projection to display it vertically? If not, is there another map source I can use that uses 2163 or can be projected to it?

If all of this is out of a question, is there another way to generate equal area grids for Mercator 900913 which will line up vertically?


texts: Read integer field as a string for list comprehension in arcpy

texts:
I have a script that is used to update a line feature class with the start and end elevations based off of a shared ID number with a point feature class using list comprehension in arcpy.  However I have an issue where the ID field from the points is a long integer and the ID field in the lines is a string field.  Can I write something in my script to read in the integer values from the point layer as strings so that the list comprehension will work?  Here is my current script that runs between matching integer field types, but not between mismatched field types:

from time import strftime
print "Start script: " + strftime("%Y-%m-%d %H:%M:%S")
import arcpy

#Transfer of Multiple Field Values between Feature Classes where there is a 1:1 Match between Field Sets
#keyValue fields must be identical field types (string, integer, etc.)

sourceFC = "Boxes"
strBoxID = str('Box_ID') #read the 'Box_ID' field as a string due to mismatching fieldtypes between boxes and lines

sourceFieldsList = [strBoxID, 'Elevation', strBoxID, 'Elevation']  # Your Source fields

# Use list comprehension to build a dictionary from a da SearchCursor
valueDict = {r[0]:(r[1:]) for r in arcpy.da.SearchCursor(sourceFC, sourceFieldsList)}

updateFC = "LinesTest"

updateFieldsList = ['FBOX', 'F_ELEV', 'TBOX', 'T_ELEV']  # your updateFC fields

with arcpy.da.UpdateCursor(updateFC, updateFieldsList) as updateRows:
    for updateRow in updateRows:
        # store the Join value of the row being updated in a keyValue variable
        keyValue = updateRow[0]
        # verify that the keyValue is in the Dictionary
        if keyValue in valueDict:
            # transfer the values stored under the keyValue from the dictionary to the updated fields.
            for n in range (1,len(sourceFieldsList)):
                updateRow[n] = valueDict[keyValue][n-1]
            updateRows.updateRow(updateRow)

del valueDict

print "Finished script: " + strftime("%Y-%m-%d %H:%M:%S")


texts:
I am a regular user of ArcGIS for Desktop (Basic License). Every now and again a drawing will have a green dot appear, see screenshot;

Sometimes this dot has its coordinates beneath it (as it does here), sometimes it will not. I am unable to select this dot, I cannot find it in any of my layers, and seem unable to remove it.

I have further not been able to identify any pattern of when such a dot appears.

What is this dot, and how can I remove it?


texts: Fixing ogr2ogr without spatialite support?

texts:
I am trying to use ogr2ogr to convert OSM data from pbf to spatialite. I used to do so with the following command:

ogr2ogr -f "SQLite" -dsco SPATIALITE=YES png-full-last.db planet_137.084,-15.689_157.615,1.255.osm.pbf

I am using Ubuntu Trusty 14.04. Recently I did an attempt to update GDAL to 2.0.1, which did not worked out. Later I attempted to completely remove GDAL using Synaptic. For this I removed also gqis, python bindings, gdal-bin, python-gdal, libgdal-dev and libgdal1h. Later I removed, clean and re-added all repositories. For some reason even when everything is removed, ogr2ogr is still available in the terminal (and I could not find a way to remove it).

Later I tried to re-install qgis with apt-get. I still get the same result, which is the following:

ERROR 6: OGR was built without libspatialite support
... sorry, creating/writing any SpatiaLite DB is unsupported

Can anyone orient me on how to fix ogr2ogr so it has again spatialite support?


texts:
I have studied about ArcSDE schema and SDE tables in Oracle, MSSQL and also PostgreSQL. But after my research I still have many questions that I haven't been able to find out.

Initially I want you to take a look at this image.
http://help.arcgis.com/en/arcgisdesktop/10.0/help/002n/pdf/sdesystables_diagram.pdf

In this diagram you can see SDE tables, but the thing is we don't know when to reconcile/post/edit also move -whatever I want- how does this mechanism work?

As you can see on the diagram, there is no relationships in these tables. In fact they are, but we don't know what kind of relationship.


texts: How to display SWPC aurora forecast data in Leaflet

texts:
The US Space Weather Prediction Center publishes the output of its Ovation forecast model on their website as JSON files: https://services.swpc.noaa.gov/json/ovation_aurora_latest.json
As far as I can tell, the file is pretty straightforward but it's not a standard GeoJSON file. I'm a beginner both with Leaflet and with JavaScript.
What would be a good way to get this data as an overlay to a Leaflet map?
Should I try to knit a function which reads the JSON, creates a new GeoJSON layer and creates a bunch of polygons or would there be a smart way to pre-convert the data somehow to some format I can just load into Leaflet?
I have touched GDAL here and there but I'm far from "fluent" with it.  I imagine, doing the conversion in the browser in JavaScript would be fairly expensive.


texts: Using multiple workspace environments ArcGIS Python

texts:
How do I combine all three workspaces?

Only 1 workspace is being read while exporting data.

import os
import arcpy

# Set environment settings
arcpy.env.workspace = "C:/eis/eis shape/UPDM.gdb/P_PipeSystem"
arcpy.env.workspace = "C:/eis/eis shape/UPDM.gdb/OtherCompanies"
arcpy.env.workspace = "C:/eis/eis shape/UPDM.gdb/NEO"

arcpy.env.overwriteOutput = True

# Set local variables
outWorkspace = "C:/data"

try:
    # Use ListFeatureClasses to generate a list of inputs
    for infc in arcpy.ListFeatureClasses():

    # Determine if the input has a defined coordinate system, can't project              it if it does not
    dsc = arcpy.Describe(infc)

    if dsc.spatialReference.Name == "Unknown":
        print ('skipped this fc due to undefined coordinate system: ' + infc)
    else:
        # Determine the new output feature class path and name
        outfc = os.path.join(outWorkspace, infc)

        # Set output coordinate system
        outCS = arcpy.SpatialReference('WGS 1984')

        # run project tool
        arcpy.Project_management(infc, outfc, outCS)

        # check messages
        print(arcpy.GetMessages())

except arcpy.ExecuteError:
print(arcpy.GetMessages(2))

except Exception as ex:
print(ex.args[0])


texts:
I have an excel file with 2 x,y coordinates as columns. One set of coordinates is the top left vertex of a rectangle, the other set is the bottom right vertex of the rectangle.

Is there anyway to plot these coordinates and have Arcmap automatically draw rectangles connecting these vertexes? I've searched on here for the answer but kept getting hung up on some of the options. Preferably a solution without using arcpy.


texts: Retrieve image statistics as numbers in Google Earth Engine

texts:
I am working with a simple piece of code for obtaining EVI index from two Landsat8 images and then calculating the difference between both.

var img_t0 = ee.Image('LANDSAT/LC8_L1T_TOA_FMASK/LC82330892017001LGN00');
var img_t1 = ee.Image('LANDSAT/LC8_L1T_TOA_FMASK/LC82330892017081LGN00');

function EVI(image) {
  var step1 = image.expression(
    '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
      'NIR': image.select('B5'),
      'RED': image.select('B4'),
      'BLUE': image.select('B2')
      });
  var step2 = step1.updateMask(image.select(['fmask']).eq(0));  // mask everything different than 'clear'
  return step2
}

var n0 = EVI(img_t0);
var n1 = EVI(img_t1);

var diff = n0.subtract(n1);

// Mean and standard deviation reducers.
var meanReducer = ee.Reducer.mean();
var sigmaReducer = ee.Reducer.stdDev();

// Use the reducer to get the mean and SD of the image.
var mean = ee.Number(diff.reduceRegion({
  reducer: meanReducer,
  bestEffort: true,
}).get('constant'));

var sigma = ee.Number(diff.reduceRegion({
  reducer: sigmaReducer,
  bestEffort: true,
}).get('constant'));

print(mean, sigma, mean+2*sigma)

From the resulting image, I need to get the mean and standard deviation for further calculations. Despite the fact that I can print the parameters (as numbers), I can not retrieve their values as such type. When I use them in an arithmetic operation, I am getting NaN as result.

I understand that the reducer property gives you a dictionary as a result, but I have tried several methods (link casting the variables with ee.Number or using brackets [ ] for getting the key value without success).

So the question is ¿How can I retrieve those values like numbers or additional variables for using them later?


texts: Simple script to batch-convert KML to SHP or G

texts: I'm not sure if this is possible given the looseness of KML files, but I have a couple thousand KML files that I would like to convert to SHP files. They are almost exclusively line or point data. I believe it is mostly possible to convert kmls via a GUI functionality in QGIS or ArcGIS but I would hopefully be able to do them all at once via command line.


texts: Automatically save python scripts in QGIS

texts:
When writing python scripts, QGIS does not save them when the user saves his main project with the ctrl+S command or with Project->Save. Besides, when leaving QGIS, the program does not warn about unsaved changes in python scripts even if it warns about unsaved changes to the project and offers to save them.

Is there any option to change this behaviour?

I should add that using the dialogue option to save the project when quitting does indeed save the project, but not the python scripts. Those are simply lost.


texts:
What I'm working on:

I'm building an app that is mostly related to geodata. I consider having a lazy filing database, so it will be filling while the users add new cities to their profiles. I'm looking for a few geo services that can provide geo services for that purpose.

What I'm looking for:

autocomplete (search city by name)
nearby search (search city by lat/long)
get full details of the city (name, country, popularity/importance etc)

What I have already checked:

Geonames (not bad solution, but data structure is not accurate, autocomplete/nearby-search works not so good)
Google Maps (it's ok, but there is not enough data of the city, so I can't even get some fields like popularity or importance)

Notes:

Mostly I consider a free solution (for a start), but it's ok to have a paid solution if the price is ok.

Edited:

According to the last suggestions, I got the point:

  if you want to have a free solution - do it yourself (prepare your own geo server)

Actually it's ok for me, but I consider to have it in the future, but for now, I'm looking for a fast solution (already hosted search engine), so I don't need to play with that on my side.

So what are the best search engines to start with (without hosting it locally)?


texts: Unable to map points on a map drawn using shapefile

texts:
I am new to this area. I have a dataset with latitude longitude data of the city of London. I would like to draw the corresponding map of the city with this points superimposed onto it. I am using python for the visualization. The shapefile I downloaded is from here:https://data.gov.uk/dataset/6cdebf5d-c69b-4480-8c9c-53ab8a816b9d/statistical-gis-boundary-files-for-london . But when I produced the maps, I noticed the there is a huge difference between the lat/lon values I have and the map generated by the shapefile. Below are the figures generated:

Below is the lat/lon values I have.

It seems to me that there must be some scaling that needs to be done to put the red points on top of the map. But I have no idea how or why or if something else is wrong.

Being a complete newbie in working with lat/lon values, I am probably making a very rookie mistake. Below is also the Python code I am using:

from shapely.geometry import MultiPoint, Point, Polygon
crs = {'init':'epsg:4326'}

geometry = [Point(xy) for xy in zip(df['lon'],df['lat'])]

state_map = gpd.read_file('../statistical-gis-boundaries-
    london/ESRI/MSOA_2004_London_High_Resolution.shp')

geo_df = gpd.GeoDataFrame(df, crs = crs, geometry = geometry)

fig,ax = plt.subplots(figsize =(12,10))
state_map.plot(ax=ax)
geo_df.plot(ax=ax,markersize=20,color='red',alpha=0.4)


texts: Add OpenLayers map on map canvas with PyQGIS

texts:
I have map canvas with two shapesfiles, and I want to add OSM map behind. My code for this is:

from PyQt5.QtGui import *
from PyQt5.QtWidgets import QHBoxLayout, QDockWidget, QWidget
from osgeo import gdal
from qgis.core import *
from qgis.gui import *

layerInput = QgsVectorLayer('path to shape file1', "file 1", "ogr")
layerInputShape = QgsVectorLayer('path to shape file2', "file 2", "ogr")

# canvas settings
canvas = QgsMapCanvas()
canvas.setCanvasColor(Qt.white)
canvas.enableAntiAliasing(True)

#defining OSM
xml = """<GDAL_WMS>
<Service name="TMS">
<ServerUrl>http://tile.openstreetmap.org/${z}/${x}/${y}.png</
ServerUrl>
</Service>
<DataWindow>
<TileLevel>18</TileLevel>
<TileCountX>1</TileCountX>
<TileCountY>1</TileCountY>
<YOrigin>top</YOrigin>
</DataWindow>
<Projection>EPSG:4326</Projection>
<BlockSizeX>256</BlockSizeX>
<BlockSizeY>256</BlockSizeY>
<BandsCount>3</BandsCount>
<Cache />
</GDAL_WMS>"""
vfn = "/vsimem/osm.xml"
gdal.FileFromMemBuffer(vfn, xml)
rasterLyr = QgsRasterLayer(vfn, "OSM")

# puting both vector layers on canvas
extent = QgsRectangle()
layers = []
QgsProject.instance().addMapLayer(layerInput)
extent.combineExtentWith(layerInput.extent())
layers.append(layerInput)
QgsProject.instance().addMapLayer(layerInputShape)
extent.combineExtentWith(layerInputShape.extent())
layers.append(layerInputShape)

#puting osm on canvas
QgsProject.instance().addMapLayer(rasterLyr)
extent.combineExtentWith(rasterLyr.extent())
layers.append(rasterLyr)

extent.scale(1.1)
canvas.setExtent(extent)
canvas.setLayers(layers)

# for legend
root = QgsProject.instance().layerTreeRoot()
bridge = QgsLayerTreeMapCanvasBridge(root, canvas)
model = QgsLayerTreeModel(root)
model.setFlag(QgsLayerTreeModel.AllowNodeReorder)
model.setFlag(QgsLayerTreeModel.AllowNodeRename)
model.setFlag(QgsLayerTreeModel.AllowNodeChangeVisibility)
model.setFlag(QgsLayerTreeModel.ShowLegend)
view = QgsLayerTreeView()
view.setModel(model)
view.show()

# embed the view and canvas into dock widget
layout = QHBoxLayout()
layout.addWidget(view)
layout.addWidget(canvas)
docked = QDockWidget()
docked.setContentsMargins(9, 9, 9, 9)
dockedWidget = QWidget()
docked.setWidget(dockedWidget)
docked.setWindowTitle('title')
dockedWidget.setLayout(layout)
docked.show()

# canvas properties
QgsProject.instance().addMapLayer(layerInput)
canvas.setExtent(layerInput.extent())
canvas.setLayers([layerInput])

canvas.refresh()
canvas.show()
app.exec_()

But on the map canvas, I just get shapefiles, without OSM map.


texts: QGIS Atlas: Filter maps based on attributes (more than one)?

texts:
I have a series of maps with the same geographic boundaries but with different set of points (from the same layer).
I would like to create an Atlas where the change is based on the attributes.

A table for example:

| ID | Party | Name    | Year |
|  1 | Rep   | Trump   | 2016 |
|  2 | Dem   | Clinton | 2016 |
|  3 | Dem   | Sanders | 2016 |
|  4 | Dem   | Obama   | 2012 |
|  5 | Dem   | Romeny  | 2012 |

I would like to autogenerate 4 maps, each shows the location of the relevant representatives by their Party + Year.


texts: Importing multiple shapefiles into single PostGIS table

texts:
I have a folder full of shapefiles that I want to import into a single PostGIS table. This table should have the following columns:

id, filename, multipolygon

I can bulk import all shapefiles using the SPIT plugin for QGIS, but each shapefile has different attributes, thus they are placed in different tables and I cannot append them into just one. I don't really want any of the attributes, I just want each polygon to be associated with the file from which it came.

I've tried two basic strategies, and neither looks very promising to me:

Import all files normally and use SQL queries to extract the information I want from each table. My SQLfu is not strong enough to do this yet. In particular, how do I get the name of the file associated with the corresponding rows?
Programmatically generate .dbf files with the associated filename as an attribute. Then import all files normally and append to a single table. Seems quite hacky.

I realize very similar questions have been posted before, but none of them seem to address this issue. I am just starting with GIS tools


texts: Placing all labels in ArcMap non-overlapped by extending arrow/line from label to point?

texts:
I was wondering if there is a way to place all labels on a map non overlapping, but by extending an arrow/line from the label to the point.

I've tried the call out box, but I couldn't find a way to enable all labels, and second extend the lines


texts: How can I use QGIS to create a metro style map

texts:
I want to create a map in the style of a London tube map (or similar metro systems) where the stations are linked together (with train lines) but the placement of the stations does not directly relate to their position in real life.

In QGIS I can create the stations (as points) and can draw lines connecting the stations. How do I now:

Move the stations into different positions (i.e. away from their real life position) ensuring that the lines connected to a given station move as the station is moved?

I foresee having a layer(s) showing the original positions of the stations and their connections plus another layer(s) showing where the stations have moved to which I can modify to make the tube map clear.


texts: tauDEM flow accumulation (D8 or D-infinity) results is NA-values for some major rivers

texts:
I am having some trouble with tauDEM that I am hoping someone can give me some advice on. I'm running tauDEM via command line through R-studio. I am working with some 30m resolution SRTM DEM data in geographic projection to do some very simple hydrology GIS:

DEM pit-filling > flow direction > flow accumulation

In doing this though, I end up with NA's (not 0's) all the way down some (but only some) major rivers.

I tried burning in one of the streams since it is happening in along a very flat meandering/braided river floodplain in one case. This did not help. I tried running D-infinity algorithms instead of D8, on both a burned/non-burned DEMs and got the same result. If I bring in the resulting flow accumulation to ArcGIS and zoom in on it, I can SEE the river clearly as a line of NA's!

I don't know what is going on. I don't think it is running out of bit depth for the number of accumulated cells since it is being stored as a 32-bit raster (I don't know how to control this either). The DEM is not very large in file size (~48mb) and the process is completed very quickly on my Windows 10 PC:
Dell Precision 7510, Intel 4-core XEON CPU @2.9GHZ (overclocks to 3.4GHZ as needed), 64GB ram, 1TB SSD, 2TB HDD.

I am hoping someone has encountered this problem and can tell me what is going on. Links to the data are just below, followed by my R code. I attached a couple plots of the resulting log(flow accumulation) including one with a missing major river plotted in blue.

link to DEM data

link to zipped shapefile

#Set working directory
setwd("enter_WD_with_data_here")

#Load some packages
require(rgdal)

require(raster)

require(sp)

# Get DEM data and a shapefile
DEM <- raster("DEM.tif") # DEM is in geographic coordinates

Rio_Majes <- shapefile("Rio_Majes.shp") # this is a shapefile I delineated of one of the missing rivers

#Pit-fill the DEM with tauDEM
system("mpiexec -n 8 pitremove -z DEM.tif -fel DEM_fill.tif")

filled_DEM=raster("DEM_fill.tif")

# D8 flow directions with tauDEM
system("mpiexec -n 8 D8Flowdir -p 30m_anp.tif -sd8 30m_sd8.tif -fel DEM_fill.tif",show.output.on.console=T,invisible=F)

p=raster("30m_anp.tif")

sd8=raster("30m_sd8.tif")

# D8 contributing area in tauDEM
system("mpiexec -n 8 AreaD8 -p 30m_anp.tif -ad8 30m_ad8.tif")

ad8=raster("30m_ad8.tif")

anp_d8=raster("30m_anp.tif")

plot(log(ad8)) #this plot is pretty chaotic, but it tells us a bit

#Focal means plot to highlight problem of missing streams with flow accumulation
avg_mat <- matrix(1, ncol=9, nrow=9)

avg_facc <- focal(ad8,w=avg_mat,fun=mean, na.rm=T, filename="avg_facc.tif") #this plot may seem dumb, but "avg_facc" is used for subsequent analyses

avg_facc=raster("avg_facc.tif")

plot(log(avg_facc)) #here is a decent plot of drainages

par(mfrow=c(1,1), new=TRUE) #this lets us add the shapefile to the plot in the right place. The plot looks good in  the console, but messes up if you "zoom" in Rstudio. Exported plots are fine.

plot(Rio_Majes, add=T, col="blue", lwd=1) #with this shapefile added we see where a river channel was/is absent

#Plot of focal average Log D8 flow accumulation

[![enter image description here][3]][3]

#Plot of focal average Log D8 flow accumulation with shapefile added


texts:
I am performing a pretty simple query in PostGIS 2.1.1 from within QGIS 2.2 for Windows (using  the DB Manager plugin). I have a table with many columns and I am trying to perform a Group By function. However, PostGIS returns an error when I do not add all the columns from my SELECT call within the GROUP BY call (see below).

ORIGINAL QUERY

SELECT city, state, city_abbr
FROM "public".footprints
WHERE city_abbr = 'nyc'
GROUP BY city

This query returns the following error:

column "footprints.state" must appear in the GROUP BY clause or be used in an aggregate      function
LINE 1: SELECT city, state, city_abbr

If I change my query to the following, it runs.

ALTERED QUERY

SELECT city, state, city_abbr
FROM "public".footprints
WHERE city_abbr = 'nyc'
GROUP BY city, state, city_abbr

I do not want to have to add each column to my GROUP BY statement. Can anyone shed any light on why this would be happening.

Thanks...

EDIT

Here is an example of the data in the table

gid   | city      | state   | city_abbr
______________________________________
0        boston      ny         bos
1        manhattan   ny         nyc
2        brooklyn    ny         nyc
3        newark      nj         nyc
4        denver      co         den
5        seattle     wa         sea


texts: QGIS Vector Grid tool differences in Apple and Windows?

texts:
Seeing a different Vector Grid tool in Apple and Windows OS. Apple is using 2.14.8 and Windows is using 2.16.3
Is this due to different versions of QGIS or is it an OS issue?


texts: In QGIS I want to transform a point shapefile from WGS84 to UTM32 (EPSG:4647). The result moves the point from somewhere in Norway to a position off the coast of Congo, Africa. Can anyone explain what's wrong?


texts:
I have a world DEM raster in QGIS in EPSG:4326 with elevation units in meters. I'm trying to get a raster of the slope of the terrain, and I'm doing so using the "Slope" tool under "Raster terrain analysis" inside of the Processing toolbox on the DEM.
At first I tried it with the default settings (z = 1), and it didn't work at all (areas either had a slope of 0 or a slope of around 89). Looking through this stackexchange, I found another similar question posted from before (Getting incorrect slope values from an ASTER DEM in QGIS?), where the solution (to my understanding) was to basically set z equal to 111120 or warp the raster to a projection that has elevation in meters already.
Setting the z equal to 111120 still hasn't worked for me, I'm assuming because the person in that question had a raster of a very specific location where setting z equal to 111120 is appropriate, and my raster is of the world and you can't do that.
I think what I need to do is to warp my raster into a projection that has elevation in meters, is that correct?
And if so, what projection should I choose and how should I do so?
The projections suggested in the answers to the previous question seemed to be specific to that specific location.
I'm also slightly suspicious that perhaps the elevation data in the DEM is somehow formatted wrong or something, as this is my first time working with DEMs (if that wasn't obvious already lol). Generally when working with raster analysis tools (Zonal Statistics, etc.) in the past, I have the option to choose which band of the raster to use, even if there's just one band, and I still see the option now to choose the elevation band when looking through the other raster tools.
However, I'm not seeing an option to select a band when using the slope tool- is this normal?


texts: Changing color according to attribute from another layer using PyQGIS?

texts:
I have township polygons on a first layer and landmark points on a second layer.

I would like change the color of each township polygons when landmarks points are missing.

I can already select one polygon according its code think to the following script :

import processing
layer = processing.getObject('township')
query = '"code" = 99 '
selection =    layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))
layer.setSelectedFeatures([k.id() for k in selection])

I can also change the color of all polygons :

symbols = layer.rendererV2().symbols()
symbol = symbols[0]
symbol.setColor(QColor("orange"))
qgis.utils.iface.mapCanvas().refresh()
qgis.utils.iface.legendInterface().refreshLayerSymbology(layer)

The attribute table of landmark layer contains the numero , geometry and township's code for each mark. Unfortanutaly some geometries are missing !

I need help to make a final script on PyQgis to change the color of only township polygons without landmarks geometry.


texts: Can an ArcGIS Online Graph be exported to a graphic?

texts: I have created a series of Graphs which are displayed within an AGS Dashboard. Is there any way to export the Graph to either pdf or a graphic format? I could of course make a screenshot, but many end-users don't know how to do this.


texts: Transforming Google Maps overlay to fit Ordnance Survey projection

texts:
I am currently in the process of creating a map with some administrative boundaries on it. The base map I am using is an ordnance survey map 1:10000. I downloaded an administrative border as a .kml and converted into an .svg file. When I tried to overlay this border on top of the ordnance survey map, it just wouldn't fit.

I slowly realised that this must be because of different projections between the ordnance survey map and the .kml which was obviously based on a Google Map.

Could someone point me in the right direction of how to warp the boundary overlay so it will fit the map?

I am using Photoshop to create this map so I apologise in advance if this is the wrong website to ask for help.

Or if someone could tell me how the ordnance survey grid compares to the Google Maps grid shape-wise (is it a trapezoid or an arc, etc.), I could try and warp it through trial and error.


texts:
I'm trying to create an app that contains a small map that should always work without internet connection for four different zoom levels.

To do this I know that I have to generate map tiles in png (or similar) to be used later on by Leaflet (or OpenLayers/Mapbox GL JS). This is no issue, because using Maperitive I can generate map tiles with given bounds and min and max zoom levels.

My problem is that I need to use a different style for roads, buildings, etc. To edit this style I used Mapbox Studio online editor and generated a JSON file with my custom styles.

What I need to do is apply this style before generating the resulting tiles, but I don't know how to do this.

Below there's a diagram for what I'm thinking about.


texts: Calculating Raster Stats inside Polygon?

texts:
I have a polygon (vector), which represents a very large area, and a raster file which also covers a large area.

Using QGIS I want to calculate what percentage of the vector polygon is covered by the raster file for certain level? Does anyone know how to do this?


texts: How to use v.net.connect & v.net.distance in QGIS graphical modeler?

texts:
I have been attempting to run the v.net.distance GRASS algorithm on vector layers in QGIS (2.18.12, on Windows) unsuccessfully. I assume that to prior to using this algorithm, you must connect each input to a network (e.g. road) layer (using v.net.connect), and take these inputs for the v.net.distance algorithm. For instance, using similar inputs mentioned in this tutorial, the workflow using my data should be:

In line with the tutorial example, the 'Electric Cable' serves as the 'streets', the 'Generator' as the 'hospitals' and the 'Poles' as the schools. The algorithm currently returns no output, and no error report is being generated under processing > history, which makes it impossible to debug. The error currently is that v.net.connect returns no output, which as I understand it is a must input for v.net.distance. I read this post, but I don't think in my case the advanced parameters are required and/or are the issue here.

My goal is to get a line layer that has the minimum paths from the 'Generator' (one feature point layer) to each pole (multiple feature point layer), via the Electric Cable network (line layer). I already achieved this through plugins (e.g. Road Graph plugin) easily, but I'd like to run in the modeler.

UPDATE* I have been able to complete this workflow in GRASS, but not in QGIS.

Here is an illustration of what the Road Graph plugin does that I'd like to get recursively done with the modeler.
Green points - Poles
Blue Line - Electric Cable
Red Star - Generator
Red Line - Shortest path


texts: Must I use "date" data type to animate data by year?

texts: I'm trying to animate the opening and closing of businesses by year but my data column is set to "number" not "date".  The animation is not working in that all the data points show up at once, rather than appear and disappear according to the year.  I tried changing the number of steps to correlate with the number of years but that hasn't worked either.  I'm wondering if the data type is the culprit.


texts:
I have a point file that has a field in it that lists all the streets at an intersection, eg) 1st St & A blvd & 2nd Ave etc

I want to be able to separate each street out individually into its own separate field.

I have some Python code I managed to field to use in field calculator that seems to work until it hits a field that doesnt have a 3rd street (and thus a third '&' symbol) so its erroring out.

How can I make it pass over these instances so it can run successfully? I think its erroring out on the 3rd passes when it comes across field with only 2 streets, eg 1st Ave & 2nd Ave

Code is below

Pre-Logic

def getField(inField):
  if inField is None:
    return None
  else:
    return inField.split(' & ')[0]

getField ( Names! )

Edit

This is basically what I want to do, a simple formula in excel seperates them, just trying to find the Python equivalent in ArcMAp

https://ibb.co/rGHVvT1


texts: How can I overlay text on a map, without a marker or textbox, using Folium?

texts:
I want to place text anywhere, as though it were part of a tile, no text box around it, and "shown on opening" (ie, show=True).

I don't see any text overlays in the documentation: https://python-visualization.github.io/folium/modules.html

Pure Leaflet appears to allow something like the "sweet static label" on display here:
http://leaflet.github.io/Leaflet.label The code for this is:

var m = L.marker([-37.785, 175.263], {draggable:true}).bindLabel('A sweet static label!', { noHide: true })
        .addTo(map)
        .showLabel();

That's better. But I would like text only.


texts: Ways and Vertices tables empty after osm2pgrouting tutorial

texts:
I've been following the osm2pgrouting tutorial, have compiled osm2pgrouting with necessary boost and expat libs installed, and built a postgresql db with support for postgis and pgrouting.

Importing the osm2pgrouting sample data appears to work in the terminal, and after import the 'classes', 'nodes', and 'types' tables look fine, but 'ways' table is empty and there are no vertices_ tables.

Has anyone encountered this issue before or have any suggestions?

Thanks,
Chris


texts: Cycle through / add name of clipping feature to output of clips in Model Builder?

texts:
I have a feature class named Study Areas with multiple polygons, named in an attribute field as Study Area A, B, C, etc.

I will be using these polygons to clip various polygon feature classes that represent different landscape attributes, like Land Cover, Soils, etc.  The end goal is an analysis of how much of each type of Land Cover acreage is within Study Area A, B, and so on.

I would like to do this through Model Builder (I don't know Python) but need to find a way that I could create outputs that have the Study Area name appended onto them so that the outputs can be easily organized (there's going to be hundreds of outputs).

So running through the model, I would need 1) a component that cycles through the Study Area feature class so that only one study area polygon is used for clipping at a time and 2) a clipped output named something like Study_Area_A_Land_Cover

Thanks


texts: Returning feature with maximum OBJECTID from Select By Attribute tool of ArcMap?

texts:
I want to get the maximum OBJECTID from a layer using "select by attribute" tool. But the result is not with the maximum OBJECTID.

When I use "select by attribute" tool to get the maximum last_edited_date, I also cannot get the latest date value.

The feature with the maximum OBJECTID and last_edited_date is newly created.

The sql statement used in the :select by attribute" is as follows:

The result of the selection:

The reason of getting the maximum OBJECTID and last_edited_date is I want to set the attribute value for the feature and after that create a buffer around the newly created feature.

I am using arcmap 10.2.


texts: GeoTiff import appears to be flat in BlenderGIS

texts: I an very new to GIS and I am trying to import a GeoTiff file into Blender to view it as a model before importing into unity. For some reason when I import the file into blender as a raw DEM it appears to be flat.


texts: ArcGIS geoprocessing DomainToTable_Management script crashing

texts:
I needed to dump all of the domains in my file geodatabase so missing values can be identified. The script runs fine if the domain name has no spaces but, unfortunately this one has many multi-word domain names. Renaming the domains isn't practical and ModelBuilder is not an option. I can list (ListDomains) all of the domains, regardless of the name and loop through everyone of them. The  Toolbox version successfully but there are too many domains for that. Here is the script.

import arcpy

domains = arcpy.da.ListDomains(pathname)

fileout = open('filename', 'w')

for domain in domains:

    name = domain.name

    print(name)

    txtline = name,'\n'

    fileout.writelines(txtline)

    arcpy.DomainToTable_management(inTable,name,outTable + name,"Code","Description")

fileout.close()


texts:
About ArcGIS Collector app, what's the difference between accuracy in the setting tab and horizontal/vertical accuracy as indicated left-bottom on the map while surveying?

For example, if I set the accuracy at 10cm with 95% in the setting menu, but I got 7.5m horizontal, 11.4m vertical accuracy when measuring a point. I don't know the different between them


texts: What's the proper way of deleting versioned data from an SDE geodatabase?

texts:
I haven't really read this anywhere.

What is the proper way of deleting versioned data from an SDE geodatabase to make sure nothing is left behind in the A & D tables?

Do I need to unversion the data before I delete it? Do I need to remove the user versions? Or does it matter?

Thanks.


texts: Exclude certain values in string field when labelling in ArcGIS Pro

texts:
I am trying to label a field on a map but exclude certain text strings from appearing. Specifically, I am working with pipes and their labels. The string values include the pipe size and material so one valaue could be (6" PVC) for example. Some of the data has 0" as the label and I am trying to remove that from the actual labeling.
Is there a way to write a python script that excludes the 0" string values from labeling but keeps all the other values?
The field name is [PIPE_IN_MA].

I decided to go with VB script and figured out what my mistake was. I only had one quoation mark when I shouldve had two. Here is the finished code:
Function FindLabel ([PIPE_IN_MA])
If [PIPE_IN_MA]="0""" Then

          FindLabel = ""
Else
FindLabel = [PIPE_IN_MA]
End If
End Function


texts: We are in the middle of creating GIS curriculum based on QGIS. When we started in June 2014 it was still 2.2, then 2.4, now 2.6.  You think it's worth returning to step 1 and begin again in 2.6?


texts: Extract longitudinal pixel values from raster.list, save to data frame

texts:
I have a list of rasters of the same location for multiple years. The change in pixel value over time represents the time series of pixel.
To further analyses, I need to extract the values over time per every pixel, and store it in data frame, where row = #pixel, column = year

Dummy data:

library(raster)

# create raster data from scratch
# create empty raster
y1<-raster(ncol = 3, nrow = 3)

values(y1)<-1:9

projection(y1)<-CRS("+init=epsg:4326")

# create and diversify the rasters
y2<-y1+10
y3<-y1+20
y4<-y1+30

# make list of rasters
y.list<-list(y1, y2,y3,y4)

# plot all rasters at once
par(mfrow = c(2,2))

for(i in 1:length(y.list)) {
  plot(y.list[[i]])
}

How the dataframe should look like:

         y1  y2  y3  y4
pixel1   1   10  20  30
pixel2
...
pixel9   9   19  29  39

How to extract unique pixel values over time, and convert individual pixel data to data frame??


texts: Create XY Event Layer using Python Script tool

texts:
I am trying to create a script toool for Making XY Event layer , But having stuck with GetParameters function while running the script with an error window of "This tool has no parameters"

Here is my code

#Import system modules
import arcpy

import os

arcpy.env.workspace = Working

arcpy.env.overwriteOutput = True

# Script arguments
Input_Table = arcpy.GetParameterAsText(0)

Output_Workspace = arcpy.GetParameterAsText(1)

# Local variables:
XY_Event_Layer = Input_Table

X Cord = "Longitude"

Y_Cord = "Latitude"

spRef = arcpy.SpatialReference(4326)

Output_Workspace = "New_Feature_Class.shp"

# Process: Make XY Event Layer
arcpy.MakeXYEventLayer_management(Input_Table, "Longitude", "Latitude", XY_Event_Layer, spRef)

# Process: Feature Class to Feature Class
arcpy.CopyFeatures_management(XY_Event_Layer, Longitude_and_Latitude_Points, "", "0", "0", "0")


texts: How to change the mouse cursor of an ArcMap add-in tool using Arcobjects?

texts:
Is it possible to change the mouse cursor of an ArcMap addin-tool using Arcobjects?  I have the cursor set in the Config.esriaddinx file, however I would like to change the cursor to reflect the underlying functionality of the tool.  The tool has select and measure functionality.  I have added a second cursor (.cur) file to my project resources in VS 2013.  I have found varying solutions, of which none have worked for me, including:

How can I access standard ArcMap cursors?

ArcGIS Add-In Custom Mouse Cursor

Set Cursor Property for ArcGIS 10 Add-in Tool?

HELP: Problems setting mouse cursor for AddIn tool

I have tried changing the build action to AddInContent, Embedded Resource, and Resource.  While using Embedded Resource, the following returns "Image format is not valid. The image file may be corrupted. Parameter name: stream"

m_customCursor = New System.Windows.Forms.Cursor(Me.GetType.Assembly.GetManifestResourceStream("MyAddin.customCursorFile.cur"))
MyBase.Cursor = m_customCursor

Using IMouseCursor.SetCursor to change to a builtin cursor seems to work, but the cursor only changes for a moment, and then returns to the cursor set in the Config file.

Dim newCursor As IMouseCursor = New MouseCursorClass
newCursor.SetCursor(2)

I have tried several of the linked esri examples for custom cursors without any success.

Is there a definitive method that can be used to dynamically change the mouse cursor for an Arcobjects addin tool?


texts: Summing attribute values for areas where multiple polgons overlap using ArcGIS Desktop?

texts:
I have one shapefile with many polygons (at least 100). Each polygon has a value. Almost all these polygons overlap with each other. I want to sum the values of the polygons in the overlappings depend on which polygons overlap in that overlapping area. The picture (simplified sketch) hopefully explains this.

But I can't find a method to achieve this in ArcMap (I use Arcmap 10.2.2). I have tried many times the Union tool and others like Merge, but I can't find a solution.


texts:
What would the best method be for determining whether an address (geocoded to a lat/long point) lies within a given shape in a GIS shapefile. Best case scenario, we are given a boolean result, that we could then apply some javascript to read yes/no results in our webapp.

We are currently using the Mapbox API for displaying the maps. I just wanted to mention that in case Mapbox has some sort of functionality in this area.


texts: Reading GeoPackage with large feature classes is very slow

texts:
I'm running a project which provides OpenStreetMap data on a global scale as GeoPackages. The data is created/transferred from a PostGIS database to GeoPackage by ogr2ogr.

Issue: Initially reading of a GeoPackage takes several minutes in QGIS, ArcGIS, or other software. For comparison, the same data can be read in under a second when stored as a File Geodatabase.

Is there a way to improve the initial read speed of feature classes stored in a GeoPackage (created by ogr2ogr, if that is relevant)?

Example data can be downloaded here: https://download.osmdata.xyz


texts: How can i join the count of points into polygon

texts:
I have two layer (point and polygon) which are kept in the same geodatabase. I want to find how many points exist in one polygon, and I googled my problem, but I could not understand how to solve my problem exactly (with this answer)

I want to perform a density analysis based on the counts of the points in each polygon, but I cannot find how many points are contained in each polygon.


texts: Multi-line labeling of a shapefile layer with pyQGis

texts:
I want to label a shp layer on two lines with pyQGis.
I managed to do it with only one field, but with two or more fields (using concat) it does not work anymore.

Here is the code I wrote:

def labelLine (lyr, fieldNom1, fieldNom2):
  # This function shows the labels of 2 fields "fieldName1" and "fieldName2" for a Line layer lyr
  label = QgsPalLayerSettings ()
  label.readFromLayer (lyr)
  label.enabled = True
  # we concatenate the two fields
  label.fieldName = "concat ('BASE:', fieldName1, '\\ n', 'L =', fieldName2, 'm')"
  label.placement = QgsPalLayerSettings.Line
  label.bufferDraw = True
  label.bufferSize = 1
  label.setDataDefinedProperty (QgsPalLayerSettings.S ize, True, True, '8', '')
  label.writeToLayer (lyr)


texts:
I am using QGIS 3.6.

By clicking on a polygon I want that an image associated to this polygon appears (as a pop-up window e.g.). The path for the pics are stored in the attribute table (string field).

Creating a generic action:

with the following action text:

open [%image%] or open [%"image"%]

doesn´t function showing the error:

With these instructions the action still doesn´t function.

What I´m doing wrong?

UPDATE: the appropriate images are store locally (no websites, no URL links). The images should appear locally and not as a website pop up in QGIS.


texts: ArcGIS GeoRSS widget source code not recognizing com.esri.ags.portal package import

texts: For some reason while attempting to import the com.esri.ags.portal.*, despite having the library imported in flex builder, it keeps on giving the error "import cannot be found." The strange thing is every other required import from the sample widget is not giving the same error. I tried multiple flex versions and refreshed and rebuilt, but it still has that issue.


texts:
I just read from GEOS 3.7 documentation that there is a new IndexedFacetDistance imported from JTS. Reading from the JTS link, I cannot tell what the new distance metric does except that it's taking advantage of R-tree index in computing something.

So, what does the distance compute functionally?

Is its result going to be the same as ST_HausdorffDistance or another existing distance measure in PostGIS, or is the distance for two given geometries going to be a different number from any existing distance?


texts: Excel table losing color in over half of the highlighted cells when exported from ArcMap 10.4 to PDF

texts:
I'm trying to add a table to from Excel into layout view. Bringing it into Arcmap it holds the color in all the highlighted cells and looks perfect. There are about 85 highlighted boxes (out of 740) to show exceedence levels of pollutants. When I export map to PDF the table loses color in about 50 of the 85 highlighted cells. I've brought it in as a JPEG but the quality is poor. I tried to insert an object and create a table but the table looks nothing like it should. It comes out very compact and cannot be read. Anyone have some tips?

Table in Arcmap

Table exported to PDF


texts: Problem with CadTools plugin in QGIS 2.2

texts:
I've just downloaded the CadTools plugin from QGIS Plugins page (can't do this through QGIS itself as I'm behind a proxy) and after installation I keep getting an error saying that it "Couldn't load plugin cadtools due an error when calling its classFactory() method".

Looking further into this there appears to be a problem with a line in the cadtools.py file relating to self.localPtah.pluginPath.  pluginPath seems to be set to look in the i18n subfolder but I don't appear to have this subfolder in my installation.

Has anybody else had a similar issue with CadTools or have I just missed a step with the installation??


texts:
I am trying to use a polygon shapefile to clip a raster, but it gives me "Clipping Error 1: Cutline polygon is invalid", and a blank raster as output.

This is similar to this question, but the answers suggested there don't seem to work. Both layers are in the same projection. I have checked and fixed the shapefile geometry, and I have tried buffering the polygons by 0m. Neither resolve the error. Is there anything else I can attempt? I was wondering about dissolving the shapefile before clipping? However, it has ~600k polygons, so I hope to avoid a step that is likely to take a long time.


texts:
I am trying to export the exact extent of a georeferenced image layer to jpg format using the exporttojpg function in arcpy. However, I am having issues setting the appropriate scale and dataframe extent prior to exporting the image. I tend to have an image export that either cuts off part of the image data or has a white boundary with no data.
My result: lots of white space surrounding image

The result I want: exact envelope export of the image

Here is the sample code I've been trying to use
    mxd = arcpy.mapping.MapDocument(mxdexport_template)
    df = arcpy.mapping.ListDataFrames(mxd,'*')[0]
    sr = arcpy.SpatialReference(4326)
    df.SpatialReference = sr
    lyrpath = os.path.join(scratch,str(auid) + '.lyr')
    arcpy.MakeRasterLayer_management(imagepath,lyrpath)
    image_lyr = arcpy.mapping.Layer(lyrpath)
    arcpy.mapping.AddLayer(df,image_lyr,'TOP')
    image_layer = arcpy.mapping.ListLayers(mxd,"",df)[0]
    image_extent = image_layer.getExtent()
    df.extent = image_extent
    arcpy.RefreshActiveView()
    arcpy.mapping.ExportToJPEG(mxd,output.jpg,df,df_export_width=5100,df_export_height=5100,world_file=True,color_mode = '24-BIT_TRUE_COLOR', jpeg_quality = 70)


texts: Options for generating travel time isochrones in QGIS 3.4: what do you use?

texts:
I'm currently using the updated OSM tools to generate isochrones in QGIS 3.4, which works well (example below).

My question is just what methods other people use and whether there is something better (e.g. live traffic information - I can dream). I have seen a variety of methods online, but I am specifically interested in methods with QGIS 3.4 - I know it can be done other ways of course but I want to make sure I'm not missing anything.


texts: How can I add to DB new created Polygon?

texts:
I'm new to LeafLet and I'm trying to understand power of LeafLet for mapping..

Question is simple, let's say I've created a new object(e.g. polygon)

map.on('draw:created', function (e) {
            var type = e.layerType,
                layer = e.layer;
drawnItems.addLayer(layer);
                });

Here, how can I pass my new object like parameter ?

I am working on android and I have to send new object to android..

One more question.. What is type of new created polygon ? How can i check it from JavaScript ?

After sending new object, I'll add to DB..


texts: Using Python to convert rasters (.asc) into rendered image GeoTIFFs (.tif) with new projection

texts:
I am trying to use the Python Console to essentially perform the same function as right clicking the layer, then going Export > Save As, and setting the following:

The raster layer is an ascii with a CRS of EPSG:28356 that I have applied styles to. I'd like to reproject it to lon/lat (EPSG:4326), and save the styling (hence needing to convert to rendered image rather than raw data). The reason I am trying to do this is because I have approximately 100 rasters (all in exactly the same location) to convert and I'd like to develop this code further to batch process all layers in my QGIS project.
I looked at a few solutions from previous questions (Solution1) (Solution2)  and pieced together their code in order to test a single layer (before making a batch):
layer = iface.activeLayer()
renderer = layer.renderer()
provider=layer.dataProvider()
pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())
width, height = layer.width(), layer.height()
extent = layer.extent()
crsDest = QgsCoordinateReferenceSystem(4326) #epsg of target crs (WGS84)
#(extra code here)
file_writer = QgsRasterFileWriter('D:/Output/outputName.tif')
file_writer.writeRaster(pipe, width, height, extent, crsDest)

I believe the problem I am having is that the extents of the original raster are in metres, while the new CRS expects degrees. I am not sure how to update the extents to suit the newly created TIFF. If I replace "crsDest" in the final line with "layer.extent()", a TIFF is successfully made, however it is in the original CRS (EPSG:28356).
In the line "#(extra code here)", I have tried two methods. For Method 1 I adopted code from another solution (Solution3). To get extents, I manually exported/reprojected the original raster with "right click > Export > Save As" (shown in the picture above) and copied the newly created extents to produce:
layer = iface.activeLayer()
renderer = layer.renderer()
provider=layer.dataProvider()
pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())
width, height = layer.width(), layer.height()
crsDest = QgsCoordinateReferenceSystem(4326)

xmin = 151.1754921309999986
ymin = -33.7468109630000015
xmax = 151.2245917170000098
ymax = -33.7199599579999969
bb = QgsRectangle(xmin, ymin, xmax, ymax)

file_writer = QgsRasterFileWriter('D:/Output/method1.tif')
file_writer.writeRaster(pipe, width, height, bb, crsDest)

This code produces a TIFF with the correct extent and CRS is produced, but is blank - it's black when viewing in an image viewing application and will load QGIS as a layer but with nothing in it. Using the info tool on the area where data should be gives this:

For Method 2, I used this solution Solution4 to produce:
import os
layer = iface.activeLayer()
renderer = layer.renderer()
provider=layer.dataProvider()
pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())
width, height = layer.width(), layer.height()
crsDest = QgsCoordinateReferenceSystem(4326)

crsSrc = QgsCoordinateReferenceSystem(layer.crs())
xform = QgsCoordinateTransform(crsSrc, crsDest, QgsProject.instance())
projected_extent = xform.transformBoundingBox(provider.extent())

file_writer = QgsRasterFileWriter('D:/Output/method2.tif')
file_writer.writeRaster(pipe, width, height, projected_extent, crsDest)

The result is the same as Method 1. A comparison between the produced TIFFs is below:

How can I make either of Method 1 or 2 work, or is there a better way to approach this problem?


texts:
I want to calculate a text field call "TIME_DAY," depending of whether or not the time falls within of four time periods it should return one of four values. My script should work, but it doesn't.

def timeday(time, am, diur, pm, noc):
    if(time >=am AND time < diur):
        return “cerp_am”
    elif(time>= diur AND time<pm):
        return “diurnal”
    elif(time>=pm AND time<noc):
        return “cerp_pm”
    else:
        return “nocturnal”

 timeday(!LMT_TIME!, !B_C_AM1!, !B_DIUR1!, !B_C_PM1!, !BEG_NOC1!)

I'm getting ERROR 000989: Python syntax error: Parsing error Syntax.

Here is what my attribute table looks like:

@Midavalo's solution worked well for the vast majority of my data, except for a few datasets, where it mistakenly places midnight values in the "diurnal" category.


texts: Creating Hillshade after selecting Raster through pGXDialog

texts:
I am trying to create an add-in button which opens a dialog box filtering only rasterdata sets. After that once I select a Raster dataset, the Hill shade operation should be applied to that Raster dataset. I have written the following code but some where I am getting an error called "Error during conversion process". Please look at the below code and let me know if any changes need to done.

Imports ESRI.ArcGIS.Geodatabase
Imports ESRI.ArcGIS.DataSourcesRaster
Imports ESRI.ArcGIS.Carto
Imports ESRI.ArcGIS.ArcMapUI
Imports ESRI.ArcGIS.GeoAnalyst
Imports ESRI.ArcGIS.esriSystem
Imports ESRI.ArcGIS.Display
Imports ESRI.ArcGIS.Geometry
Imports ESRI.ArcGIS.Catalog

Public Class AddHillshade
Inherits ESRI.ArcGIS.Desktop.AddIns.Button

Public Sub New()

End Sub

Protected Overrides Sub OnClick()
    '
    '  TODO: Sample code showing how to access button host
    '
    My.ArcMap.Application.CurrentTool = Nothing
    Dim pMxDoc As IMxDocument = My.ArcMap.Document
    Dim pMap As IMap = pMxDoc.FocusMap
    Dim bSelectedObject As Boolean
    Dim pGxDialog As ESRI.ArcGIS.CatalogUI.IGxDialog = New ESRI.ArcGIS.CatalogUI.GxDialog
    'only show raster datasets
    Dim pRasterFilter As ESRI.ArcGIS.Catalog.IGxObjectFilter = New ESRI.ArcGIS.Catalog.GxFilterRasterDatasets
    'Layers and document declaration
    Dim pGXSelect As IEnumGxObject = Nothing
    Dim pGxObject As IGxObject = Nothing
    Dim pGxDataset As IGxRasterDataset
    Dim pRasterLayer As IRasterLayer
    Dim pDataset As IDataset

    Try
        With pGxDialog
            .AllowMultiSelect = False
            .StartingLocation = "C:\Users\rxm160030\Downloads\Lab12\Data"
            .Title = "Select Raster datasets to Add to Map"
            .ObjectFilter = pRasterFilter

        End With
        bSelectedObject = pGxDialog.DoModalOpen(My.ArcMap.Application.hWnd, pGXSelect)
        'If no file selected, exit
        If bSelectedObject = Nothing Then Exit Sub
        pGXSelect.Reset()
        'Get first dataset
        pGxDataset = pGXSelect.Next
        'no dataset selected
        If pGxDataset Is Nothing Then Exit Sub
        'iterate through remaining datasets
        Do Until pGxDataset Is Nothing
            pRasterLayer = New RasterLayer
            Dim pInputRaster As IRaster = pRasterLayer.Raster

            ' Do something with the featureLayer
            Dim choice1 As Double = InputBox("Enter Azimuth")
            Dim choice2 As Double = InputBox("Enter Altitude")

            Dim pSurfaceOp As ISurfaceOp = New RasterSurfaceOp
            Dim pEnv As IRasterAnalysisEnvironment = pSurfaceOp
            Dim pWSF As IWorkspaceFactory = New RasterWorkspaceFactory
            Dim pWS As IWorkspace = pWSF.OpenFromFile("C:\Users\rxm160030\Downloads\Lab12\Data", 0)
            pEnv.OutWorkspace = pWS
            Dim pOutput As IGeoDataset = pSurfaceOp.HillShade(pInputRaster, choice1, choice2, True, 1)

            Dim pOutLayer As IRasterLayer = New RasterLayer
            pOutLayer.CreateFromRaster(pOutput)
            pOutLayer.Name = "HillShade"
            pMap.AddLayer(pOutLayer)
            ' Go to next dataset
            pGxDataset = pGXSelect.Next
        Loop
    Catch ex As Exception
        MsgBox("BtnAdd_Click Exception: " & ex.Message)
    Finally
        ' Release ArcObjects objects
        pGXSelect = Nothing
        pRasterFilter = Nothing
        pGxDialog = Nothing
        pGxDataset = Nothing
        pRasterLayer = Nothing
        ' Call garbage collection
        GC.Collect()
        GC.WaitForPendingFinalizers()
    End Try
End Sub
End Class

After clicking the add in button I am getting an option to select the Raster Dataset and input the Azimuth and Elevation values. But Later on I am getting an error.


texts: I'm trying to use the "Locate Features along Routes." I have a simple polyline and I have an excel sheet of To and From miles.  The excel sheet has no X&Y coordinates.  The excel sheet has the information on where the breaks are based on miles.  I want the geoprocessing to take care of polygons along the route and then add in the excel file.  How do I do this?


texts: Installing pre-requisites for forest classification using SCP and QGIS?

texts:
I am using QGIS 2.8.6 and I want to use SCP for forest classification from Landsat 8 oli. I  have already installed the SPP plugin.

What else do I need to install for operation of the software?


texts: Disable Mosaic Overview Information When Using Identify Tool

texts:
I have a created mosaic datasets using elevation data. When I add this data to ArcMap and use the identify tool it brings up information about both the footprint and image.

At full extent the identify information under image contains an elevation value and also the name of the overview raster within that location. When I zoom in and identify a location it will bring up the elevation value and the name of the tile with the associated attribute information.

Is there a specific function or property that I need to configure so that the identify tool only provides the elevation value?


texts: Automatisation of isochrone map creation in QGIS 3

texts:
I need to measure the area of isochrones for a big number of points, for different time thresholds and for different means of transport.

I've been using the plugin isochrones in OSM tools with QGIS 3.4.4 but it takes a long time if there are many points. It also crashes sometimes. Besides, it needs to change the threshold manually at each step.

Is there any way to automate this processing and to avoid crashing?


texts:
I need to make a query to measure the distance from a given point to all points in a table and return the five closest to that point. The problem is I can't convert the distance to kilometers.

Here is the query:

SELECT gid,name , ST_Distance(ST_Transform(g1.geom,2163),ST_Transform(ST_GeomFromText('POINT(37.791458999999996 26.698004299999997)',2100),2163)) as km
    FROM samos As g1
    WHERE ST_DWithin(ST_Transform(g1.geom,2163),ST_Transform(ST_GeomFromText('POINT(37.791458999999996 26.698004299999997)',2100),2163), 50)
    ORDER BY ST_Distance(ST_Transform(g1.geom,2163),ST_Transform(ST_GeomFromText('POINT(37.791458999999996 26.698004299999997)',2100),2163))
    LIMIT 5;

The coords is from Greece the table is in 2100 Greece WGS geom(Point,2100)
Version is 9.3 of postgresql
the out put of cords in table samos
coordinates


texts: Weird-looking France after OSM import with osm2pgsql

texts:
After importing france-latest.osm from Geofabrik, using osm2pgsql, it seems that the polygon attributed to Metropolitan France (osm_id = 1403916) is a bit odd:

I'm far to be a geography expert, but I think that the shape is a bit off here.

Is there something I missed?

Thanks


texts: How to only draw an inside or outside stroke around a polygon in GeoServer

texts:
If I add a stroke to a polygon SLD style, I get something that extends over both sides of the polygon boundary like this:

I want to highlight certain polygons, some without drawing anything over their area, others without drawing anything over their neighbors. So the stroke should only be on the inside OR outside instead like in these:

How can I achieve this in GeoServer SLD in an easy, fast and hopefully non-hacky way? I need to be able to specify the width in "drawing units", not my geometry's coordinate system units so the suggestions from the similar Line stroke offset in Geoserver do not apply.


texts:
I want to know if is possible to change spacing of a normal (not raster) legend in Geoserver:

I need to add more space between the squares and the text.

There are dx,dy,mx,my parameters but they only work for Raster Layers


texts: How can I remove the gap between two polygon boundaries in QGIS?

texts:
I have a couple of messy polygon boundaries (example data) that should share a boundary but don't:

I've had a look at some suggested tools for a similar issue using ArcGIS in a separate question but need some sort of method for QGIS (or another FOSS4G solution).

Some possibilities:

v.clean with 'snap' only seems to work on lines (so perhaps a polys-to-lines + v.clean/snap + lines-to-polys might work?);
a combination of geoprocessing steps, perhaps buffering each polygon (filling the gap) then intersecting them?


texts: Access individual bands and use them in map algebra

texts:
I have a list of multiband rasters, consisting of: Red=Band 1, NIR=Band 2, SWIR1=Band 3  and SWIR2= Band 4. What I want is to access only the red and NIR bands and save it to new file. Then, I'll be using them my Map Algebra calculations below.

(Note that accessing and saving the Red and NIR bands portion is still missing from my code)

    import arcpy
    from arcpy import env
    from arcpy.sa import *
    env.overwriteOutput = True

    #Set the current workspace, assuming that all red and nir bands are already
    #extracted and saved as individual bands
    env.workspace = (r"C:\thesis\hansenwipfiles\first and last v1.0")

    #Raster List (Assuming the rasters here are already single bands)
    RedbandList = ["red20N_120E.tif", "red10N_120E.tif", "red20N_110E.tif",
                   "red10N_110E.tif"]

    NIRbandList = ["NIR20N_120E.tif", "NIR10N_120E.tif", "NIR20N_110E.tif",
                   "NIR10N_110E.tif"]

    for range in (0,4):
        ndvi = Raster(NIRbandList[f])-Raster(RedbandList[f])/\
                  Raster(NIRbandList[f])+Raster(RedbandList[f])
        ndviFloat = float(ndvi)
        ouputName = ndviFloat
        ndviFloat.save(ouputName)

    print "Finish!"

I have not tested the code if it works, though no errors  found after running the "Test Module".


texts: Setting Google street basemap language in QGIS?

texts:
I use google street map as my project's basemap in QGIS: Using Tile Server: https://mt1.google.com/vt/lyrs=m&x={x}&y={y}&z={z}.

As my computer default language is Chinese, when the basemap is loading, all the place names were present in Chinese, how can I change it to English?
I have already set the QGIS language to English (Option->General->Override system locale->English), but it does not work in google street basemap.


texts: Drawing distances from one point to another using QGIS?

texts:
I want to be able to generate a distance line between two points on the map, which also has the distance as a label (similar to the second image seen on Google maps or in the first image from the red point to Aberdeen).I have used the measure line tool and it has provided myself with a distance, however I am yet to figure out how to produce a line with the classified distance labeled on it. does anyone know how to solve this issue. I am using QGIS versions 3.2

I want to record the distance from the red point to the location of Aberdeen so that when I go into {print composer} I can make a map which has a line running between the two distance with its distance e.g 100km recorded I don't have the distance recorded in my {attribute table} would I need to download that data then once collected, I would be able to draw the line?


texts:
i am using the arcGIS api to make a plugin for arcFM, when i try to run this code

 Type t = Type.GetTypeFromProgID("esriFramework.AppRef");
                System.Object obj = Activator.CreateInstance(t);
                pApp = obj as IApplication;

i get

System.Runtime.InteropServices.COMException(0x8000FFFF): Creating an instance of the component with CLSID {Appref CLSID HERE} from the IClassFactory faileddue to the following error: 8000ffff


texts: Seeking Septic Tank Location Data for Southeastern United States?

texts:
I am compiling a database of septic tank locations to make some rough estimates of the magnitude of return flow to the groundwater system for the Floridan aquifer system. Ideally I would have a lat/long for each tank so that I can identify, at the census block level, the average number of persons per household to compute the number of gallons of flow-through per day. A dataset of this detail has been developed for the state of Florida and possibly for Georgia as well. Records for South Carolina and especially Alabama are not as robust and I would need to resort to some county-wide tally. I have only managed to find state-wide counts on uscensus.gov through 1990. Block-group level data seem to be available through American Fact Finder 2 from the 2010 census, but only for American territories such as Guam or the Virgin Islands (Topics -> Housing -> Physical Characteristics -> Sewage Disposal).

Can anyone point me toward any datasets which are at least at the county level?

Depending on the quality I may elect to use such a dataset for the entire estimation process to maintain continuity between states.


texts: How to use the method makeDifference in the console?

texts:
I would like to do a "difference" between 2 features in PyQGIS.

I don't understand how to use the method makeDifference() as shown on the documentation here.

It says the parameter to pass the function is:

QgsGeometry makeDifference (const QgsGeometry &other)

I have tried the following code but don't see any change in the feature.

layers = iface.legendInterface().layers()
layer = layers[0]

# I select the feature that will erase all the others
features = layer.selectedFeatures()

#Get the other features
layer.invertSelection()
inverse = layer.selectedFeatures()

#Erase all the other features with my selected feature
for f in inverse:
    try:
        layer.dataProvider().changeGeometryValues({f.id():
        features[0].geometry().makeDifference(f.geometry())})
    except:
        pass


texts: How do I automatically fill a table column?

texts: I would like to solve a basic problem. I need to automatically fill a table in QGIS. Is it possible to fill a table according to this: if column "SPECIES" contains string "DOG" → write to column "CODE" integer 1?


texts: Get information from a tool in ArcObjects

texts:
I have the following situation in an ArcMap add-in that I'm building:

This add-in contains a form where a user needs to click on a button, and then be able to draw a polyline on the map (the form is non-modal). The polyline is needed in the form for further processing.

I have now implemented this using a tool that inherits from ESRI.ArcGIS.Desktop.AddIns.Tool, and uses a RubberLineClass in the OnMouseDown event. The tool is activated in the button's OnClick event. This all works, and the user can draw the polyline.

The problem is getting the polyline from the tool. I was thinking to have the tool add the polyline as a graphic element to the map, and in the form listen to the ItemAdded event on IActiveView. Unfortunately, this event is only raised when maplayers are added. It does work with elements on the PageLayout, but if I add the element to that, the polyline's coordinates are altered. Also, this feels a bit clumsy.

So, the question is: Is there a better way to get information from a tool? Or is there a way to draw a polyline on the map without using a tool?

Maybe not relevant, but here's the code anyway:

// Activate the tool
var documentBars = ((IDocument)ArcMap.Document).CommandBars;
var cmdItem = documentBars.Find(new UIDClass { Value = ThisAddIn.IDs.PolylineTool });
ArcMap.Application.CurrentTool = cmdItem;

// The tool's OnMouseDown event
protected override void OnMouseDown(MouseEventArgs arg)
{
    var rubberBand = new RubberLineClass();
    var geometry = rubberBand.TrackNew(ArcMap.Document.ActiveView.ScreenDisplay, new SimpleLineSymbolClass { Color = new RgbColorClass { Red = 255 } });
    geometry.SpatialReference = ArcMap.Document.FocusMap.SpatialReference;

    if (!geometry.IsEmpty)
    {
        var container = (IGraphicsContainer)ArcMap.Document.FocusMap;
        container.AddElement(new LineElementClass { Geometry = geometry }, 0);
    }
}


texts: How to find the polygon inside which a point lies?

texts:
I have a layer with polygonal features. Each feature has attributes and values. I also have a list of co-ordinates and I would like to know which feature (or polygon) the co-ordinates lie in.

Could someone please guide me on how to go about this? Is there a function in the API that can help me achieve my goal or should I use some computational geometry algorithm to do it myself? I know how to do the latter but it would save me some time if there was a built in function already.

Thanks.


texts: Using Linux, is there something out there to quickly view a dbase file? Sometimes there are look-up tables(.dbf) to accompany a shapefile. Libre Office Calc can't read it. Maybe a QGIS plugin that can read and or join?


texts: Create Array From Selected Features Based on the Order in Which They Were Selected

texts:
I need to create an array of features from a selection based on the order in which they were selected.

As in: I select feature 5,2,3,4. I need to have an array so I can apply attributes (with python) to to these features based on order in which they were selected. This is too apply part numbers to pre-mapped survey plans. To add a part to each is time consuming when their are 50 parts. This would allow me to select them in the order of part number, then apply the correct part with a button I have already created with python.

I've been trying for a long time to figure this out and I have run out of solutions.

I am pretty well versed in Python and with that ArcPy. Seem's like this should be easier then it is.

Thanks for any help.


texts: Can I plot points on a map using JSON data from a Node.js API?

texts: I am new to the arcGIS JavaScript API.  I have a Node API that is providing geographical information in JSON.  I have google'd for examples and searched for solutions on this web site.  Can anyone point me towards an example where JSON data is used to plot points on a map?


texts:
I am working on a watershed analysis project. The goal of the project is to create a model that will simulate changes the effect that different land uses will have on the watershed as a whole.

Using python, I need to create a user friendly model that will allow a non-tech savvy person to accomplish the following:

"Paint"/select a group of cells in a raster land use layer and then change to their desired land use. For example, the user may wish to see what effect changing a certain area from developed industrial to forest would have on the watershed.

I have tried several options and have not had much luck with any of them.

Clipping a raster based on a polygon that the user would draw around the area they want to change. The problem with this approach is that the raster layer that is clipped is not easily merged back in with the original raster layer.

Arcscan - Raster painting toolbar: this tool seems as though it should work, however, it is limited in the way it will let you select cells. Rather than letting you select individual cells to then change, it will only let you select cells that are connected. This is easy enough to do within the attribute table of a raster layer.

Does ArcMap support the ability to select individual cells within a raster layer to then change those selected cells?

Thank you,

Joe Martin


texts:
I try to export my map to picture with the same method like this example of OpenLayers 3.

The problem is whith few baselayers, the function doesn't work and firebug say :

What can I dos for resolve this problem ?


texts: QGIS crashes when intercepting QgsVectorLayer.featureAdded event and user decides to Toggle edit mode off

texts:
I am quite new to QGIS and Python plugins.
I have the following code in my plugin
def feature_added(fid):

    # GET INSTANCE OF NEWLY ADDED FEATURE
    feature = my_layer.getFeature(fid)

    # SET SOME DEFAULTS TO NEW FEATURE
    feature['Field1'] = 666
    feature['Field2'] = 777

    # UPDATE THE FEATURE IN THE LAYER
    my_layer.updateFeature(feature)

my_layer.featureAdded.connect(feature_added)

The above code executes fine when I use the QGIS edit/add polygon feature buttons:

User successfully draws a polygon using left click for each coordinate, and finishes it off by using right click.
The user then decides to back out by clicking the Pencil icon and hits "Discard" when asked to save changes.
Here is where QGIS crashes.
I also noticed that if I removed the the code that updates "Field1" and "Field2" everything works fine, it manages to exit editing mode without crashing.
What am I doing wrong? Seem like I am blocking an internal QGIS method from ending?
Oh I am using QGIS version 3.12.3-București


texts: Erase polygon from line to get line using ArcGIS for Desktop?

texts:
How can I erase a polygon from a line to get the erased line?

Conceptual image below.


texts: How to set GCP in local coordinates system to images to be processed with ODM

texts:
I have a drone survey of a local area of 200x200 meters. I have about 45 photos.
6 photogrammetric control points were placed on the ground and measured with centimeter precision from topographic instruments.
The local coordinates system is not georeferenced, nor are there any intentions to georeference it.  So I need to adjust and reference the photogrammetric model to the local coordinates system.
I am running a dockerized version of OpenDroneMap, to which I must pass a list of references between the coordinates of the control points in the local system and in the image system, for each image in which each point are viewed.
I ran a node and a web interface only to reference the images but it seems that I can only reference them in geographic coordinates.
How can I set GCPs with local coordinates, to images to be processed with OpenDroneMap?


texts:
I want to make plan of building apartments. To make it I have to know what layers I need have and what attributes they must have (for example,communications).
Maybe somebody knows is there something like it exists? Or what can help me to do it.


texts: gdal_translate 2 bands into JPG RG image

texts:
I have 2 bands of data (specifically, u and v wind vectors from GFS in GRIB format) and I'd like to transform it to JPG image into RG channels, B channel empty. How can I use gdal_translate to transform 2 bands? gdal_translate -scale -128 127 0 255 2020051500.grib 2020051500.jpg throws "ERROR 6: JPEG driver doesn't support 2 bands.  Must be 1 (grey), 3 (RGB) or 4 bands.".

I can duplicate a band with gdal_translate -b 1 -b 2 -b 2 -scale -128 127 0 255 2020051500.grib 2020051500.jpg, but can I avoid it?


texts: Thiessen Polygon delineation within a Feature

texts:
I am trying to create Thiessen (Voronoi) polygons, based on points within a defined shapefile extent.
The output should be inclusive to the feature.  The output I receive now is the rectangular area around the selected points.

Steps I am using:
1. select points of interest (inside and outside spatial extent)
2. set primary display field to FID (aka:ObjectID)
3. Environment… > Set 'Analysis Mask' to shapefile with feature (in this case, a watershed)

Now here is where the problem starts.
I have a set of instructions for arc 9.3, but I am now operating on ArcGIS 10 (Spatial Analyst installed & activated).
Also, just to clarify, I did scour the inter-webs for some number of hours- Esri support, gis.SE, Professor Google...

########

9.3.1 version:
Spatial Analyst >> Distance >> Allocation
- set 'Assign To' as point shapefile of selected data
- set cell size

10.0 methods tried:
Analysis >> Create Thiessen Polygons
-- Many variations of selected features
-- File locations, both inside and out of file geodatabase

I am open to any solution, using Arc 10.0 and/or python would be ideal.  The first image is an example.  The output does not need to be exact.

wiki: http://en.wikipedia.org/wiki/Voronoi_diagram


texts: Where to find POI data for places of worship (e.g. churches)?

texts: I am looking for free U.S. state-wide and/or county-wide shapefiles that include geographic locations for places of worship.  Examples include churches, temples, synagogues, and mosques.  In other words, I am looking for POI (points of interest) data, nothing more.  The attribute table fields would include something like FID, Shape, and Name of Place of Worship.  I am specifically looking for data in the U.S. state of Pennsylvania, but I eventually would also like the same location information for other states.


texts: Keep or increment field value based on value of another field in arcpy

texts:
My requirement is to create auto increment numbers in the field2  based on the values in field1. But if field1 having same for more than one row, Increment number also remains same.

ex:

field1         |  field2
479049.620600  |  1
479049.620600  |  1
478611.065000  |  2
478732.117500  |  3
478890.969400  |  4
478890.969400  |  4
478890.969400  |  4
479046.462900  |  5
478576.747800  |  6

I have searched and got some code but it does not meet my requirement.

import arcpy
from arcpy import env

path = r"C:\SOFTWARE\Shapefiles\shp\Points.shp"
arcpy.env = path

rows = arcpy.arcpy.UpdateCursor(path)
myList = []

i = 0

for row in rows:
    if i == 0:
        value = row.X
        myList.append(value)
    i += 1
    print value, i
    if row.ID != value:
        value = row.X
        myList.append(value)
    rows.updateRow(row)
del rows


texts: How to fit polylines on base map layer using QGIS?

texts:
I have transformed a pdf into polyline and i want it to fit into the base map of OpenStreetMap.

How can I rotate and reduce the polylines layer?

With the image will be easier to know than I ask.


texts: Calculating length of common boundary lines of two polygon in QGIS?

texts:
I have a map consisting various district. How can I estimate the length of the common boundary shared by any two corresponding districts/polygons in QGIS?

For example if a polygon "A" shares a common boundary with polygons, "B", "C" and "D". Then what would be the length of the boundary between polygon "A-B", "A-C" and "A-D"?

I have calculated the neighboring polygon using the following in field calculator

aggregate(
  layer:='file_districts',
  aggregate:=' concatenate',
  expression:= DIST_ID,
  concatenator:=', ',
  filter:=touches($geometry, geometry(@parent))
)

Now i want to estimate the length of each of the neighbor polygon.


texts: qgis-grass-plugin7 problem updating 2.18.12 ->2.18.13 (win10-linux)

texts:
after updating qgis to  2.18.13 version (and GRASS from 7.2.1 to 7.2.2) via OSGeo4W advanced setup on win10 (and via repository on linux mint 18.1) the qgis-grass-plugin7 has disappeared from the list of available plugins.

This is the processing->options->programs:

How can i fix this problem?


texts: Shapefile is not shown on a raster image although having same CRS

texts:
I have a shapefile which was digitized from a jpg map. the only info i knew from that jpg, the map was created in lambert conformal conic(north).
so, i have digitize and create my shapefile with the same CRS.
now i want to check my shapefile by overlapping on a landsat image which is in EPSG: 32645 WGS 84/ UTM ZONE 45.
i have changed my shapefile´s CRS into EPSG: 32645 WGS 84/ UTM ZONE 45 but still the shapefile is not shown on the exact location.
what can I do now?


texts:
I'm trying to create a simple application using QGIS 3 SDK. I have used osgeo4w-setup-x86_64.exe for installing QGIS, and I have also installed qgis-dev and qgis-ltr-dev and Qt5. I have set up my .pro file like this:
QT       += core gui

greaterThan(QT_MAJOR_VERSION, 4): QT += widgets

CONFIG += c++11

INCLUDEPATH += "C:\OSGeo4W64\apps\qgis-ltr-dev\include" \
                "C:\OSGeo4W64\apps\qgis-ltr-dev\bin"

LIBS += -L"C:\OSGeo4W64\apps\qgis-ltr-dev\lib" -lqgis_core -lqgis_gui -lqgis_app

SOURCES += \
    main.cpp

and the main.cpp is like this:
#include <qgsapplication.h>

#include <QDebug>
int main(int argc, char *argv[])
{
    qDebug() << "*** QGS Application Instance Created *** ";

    QgsApplication app(argc, argv, true);
    int exitCode= app.exec();

    return exitCode;
}

I run my code in this enviroment:
Windows 10
Qt 5.15.1(MSVC 2019)
Qt Creator 4.13.1
QGIS 3.x.x

This code is compiled correctly. But, the problem is when I run it, the program crash even before entering the main function. When I comment the line QgsApplication app(argc, argv, true);, the crash does not happen. I believe there is some problem with linking the QGIS to my application.
How can I fix this?
Have I done something wrong?
Haven't I considered something in setting up the environment?


texts: Proper transformation file to go from NAD83 CSRS UTM 17 to NAD83 UTM 17

texts: Looking to determine which transformation file to select when using the Project tool in ArcMap 10.2.2 to go from NAD83 CSRS UTM 17 and NAD 83 UTM 17.  The data I am working with is from Geogratis and will be projected to suit Ontario standards (NAD83 UTM 17).


texts:
Why is ArcMap 10.5.1 search window not showing properly?

Out of ca 50 times opening ArcMap, only twice or so it looked like in previous versions, i.e. there was possibility to work with that - e.g. enter a name of function and search for it...

I am using Windows 7 Professional.

I have
-renamed ESRI folder, also in registry, not helped!
-deleted Normal.mxt, no change!
-chosen "repair" in instalation file window, also not helped!


texts: How to collect shpfiles from "Iterate Feature Classes" into SINGLE input to "Intersect" tool?

texts:
Need to collect output files (Iterate Feature Classes) into single input to INTERSECT tool.

Collect values isnt working (wrong data type)!

Only precondition is available (see pic):

http://img407.imageshack.us/img407/2945/54132406.png

Any ideas?


texts: MSSQL Spatial Query (View) is Slow in ArcMap

texts:
I created the following view in Microsoft SQL Server 2008 R2:

CREATE VIEW test AS
SELECT
    p.id,
    p.shape
FROM
    p INNER JOIN r ON p.shape.STIntersects(r.shape) = 1
WHERE
    r.Name = 'Test'

The view runs very quickly within MSSQL if the spatial index is enabled on P.  The view runs very slowly within MSSQL if I disable the spatial index.

If I add the view to ArcMap 10.3.1 (via Catalog) it is very slow.  It doesn't matter if the spatial index is enabled or disabled - the view still runs very slowly within ArcMap.

How can I improve the view's performance within ArcMap?

Update:

The SQL actually being called by ArcMap seems to take the following form:

select
    SHAPE
from
    (select id, SHAPE from test) a
where
    SHAPE.Filter(geometry::STGeomFromWKB(@P1,@P2)) = 1

I think that the WHERE clause translates to SHAPE.Filter(SHAPE) = 1.  When I run this query within MSSQL it returns the results but now it's as slow as ArcMap.


texts: How can I properly clip a shapefile to get the section I want using QGIS?

texts:
I've tried the other "how to clip" instructions here, but consistently get an empty shapefile from the clip.

QGIS Version: 2.4.0-Chugiak

OS: MacOS 10.10

Input File: San Diego (South Map) 1999 Topo (2' intervals) (1 GB Download Link)

I've assigned the same CRS from the San Diego Vector Layer to the polygon clipping mask.

Here is a screen capture of my process:

QGIS Vector Clipping Path Producing Empty File

Hopefully someone can point out where I'm going wrong with this process.


texts: Where to download QGIS win7 x64 installer?

texts:
When I try to print big map (1200mm x 1200mm x 300dpi) on pretty strong computer (i7 975 12GB RAM) withj Windows 7 I obtain error Error - 'Image too big. Creation of image with bla bla pixels failed' - same on QGIS 1.7.4 and 1.8.0 versions.

Can anybody link x64 build of QGIS please - I'd like to check printing on that version? I don't wanna install programming tools and see same error after...


texts: Import CSV to QGIS 3 where the geometry is in GeoJSON format

texts:
I have a CSV file with different attributes for each polygon. One of the attributes (column name) is GeoJSON which hold the geometry. For example, in the first polygon the geojson attribute is:

{"type":"MultiPolygon","coordinates":[[[[35.1919984817505,32.58075753670527],[35.19470214843751,32.58319843158036],[35.19618272781373,32.58186046725293],[35.19124746322633,32.57771084505097],[35.19030332565308,32.57938337350581],[35.1919984817505,32.58075753670527]]]]}

I was able to import the CSV file only as a table to QGIS 3. Is there a way to tell QGIS (or possibly via pg admin 4), that the attribute GeoJSON holds the geometry but in a GeoJSON format?


texts: ArcPy Raster Calculator error " 'str' object has no attribute 'save' "

texts:
I converted a model into python script and tried to run it. It gave me an error when it gets to the raster calculator process.The message says "str" object has no attribute "save". See picture below.

My code here :

arcpy.gp.RasterCalculator_sa("\"%Reclass_Aspe5%\" + \"%Reclass_Slop7%\"", Output)

The expression: "%Reclass_Aspe5%" + "%Reclass_Slop7%"

The code it directly from ModelBuilder. Every other processes work out fine except raster calculator. Does anyone know what the problem is?

Error snippet:

  ExecuteError: ERROR 000539: Error running expression: rcexec()
  Traceback (most recent call last):   File "", line 1, in
     File "", line 8, in rcexec AttributeError: 'str'
  object has no attribute 'save'

  Failed to execute (RasterCalculator)


texts: Geopandas: Write layer back into GeoDataBase

texts:
I use geopandas (which uses fiona) to read a GeoDataBase layer (using driver="OpenFileGDB") into a geodataframe for easy manipulation.

read_shp = gpd.read_file(GDB_file, layer=layer_choice)

That works well. I then rearrange the columns and do some data manipulation until I am happy with the data.

BUT I have no idea or if it is possible to write that geodataframe back into the GeoDataBase?

I tried:
read_shp.to_file(GDB_file, layer=layer_choice)
But that just creates a new shapefile in the wrong place.

Does anyone have a solution to this problem?


texts: Converting selected polygons to lines using PyQGIS

texts:
In PyQGIS 3, how can I convert selected polygons with n sides, into another layer with n lines?

import processing

mapcanvas = iface.mapCanvas()
layer = mapcanvas.currentLayer()
features = layer.getFeatures()

line_path = "C:/Users/Me/Desktop/pyqgis/line.shp"

#### I need something like this processing method, but instead of a polygon in the input, i need the selected feature.
#processing.run("qgis:polygonstolines", {'INPUT': polygon_path,'OUTPUT': line_path}


texts: QGIS data defined properties 'offset' syntax?

texts:
So I have a property field called 'offset', type=integer, range=0-359.  The intent is to draw an svg symbol that it offset 5 units radially away from the feature coordinate at the angle (in degrees) specified in the 'offset' field.  Basically if 'offset' is 45 and the feature coordinate is 0,0 then the symbol will be drawn at x=5*acos(45) y=5*asin(45) (sorry if the trig is wrong for the moment, you get the idea).  If 'offset' is 90 and the feature is at 0,0 then the symbol will be at x=5 y=0.

SO: the question is how to write this expression in the Data Defined Properties of the marker in the style dialog box?

(After hitting Data Defined Properties in the SVG marker dialog, check the 'Offset' box and no others...)

Ignoring the formulas for a minute, to get the syntax straight, I'm trying these on an object with 'offset' value = 5.

I've tried:

"offset","offset" -> not a valid expression

'"offset","offset"' -> valid expression, but no actual offset

concat(tostring("offset"),',',tostring("offset")) -> valid expression, but no actual offset

concat(tostring(offset),',',tostring(offset)) -> valid expression, but no actual offset

Several of the attempts give 'Expression is invalid' under 'Output preview' with the following, or similar, under 'more info':

Parser Error:
syntax error, unexpected COMMA, expecting $end

Eval Error:
No root node! Parsing failed?

As a reality check, the following expressions DO work as expected:

'5,5'

concat(tostring(5),',',tostring(5))

But the following gives the 'Expression is invalid' error:
5,5

So, apparently it does want a string?  Can anyone provide the correct syntax?  Thanks!


texts: ArcMap 'go back to previous extent' doesn't

texts:
I often use the "go back to previous extent" i.e. the back button when navigating my models in ArcMap.  Alas for some reason it seems to have fixed on a certain location as my "previous extent", no matter how many places I visit in the interim.

Is this a bug or is there a bizarre setting somewhere I've triggered?


texts: Merging two line segments in ArcGIS Desktop?

texts:
I am using ArcMap 10.2 and merging several line segments with the same ID. Most of the time this works fine and I end up with one segment. (about 95% of the time)

However, occasionally the segments will merge but I end up with two separate "parts". The "parts" I am speaking of are what you can see when you "edit vertices" and look to the "Edit Sketch Properties". It says "Part 0" ....lists all vertices..."Part 1"...lists all vertices

The end goal is building "M" values for these merged segments and I haven't found a workaround other than splitting the merged segment before the second "part" starts and extend/tracing over the original part that wouldn't merge properly. This is becoming very tedious with larger segments.