texts: Edits to a spatial dataframe dont carry over when exporting to a featureclass in ArcPy, why?

texts:
I created a spatial dataframe from a featureclass. Then I edited the spatial dataframe, and exported it back to a featureclass. I was expecting the attribute table of the featureclass to keep the results of the spatial dataframe but that wasn't the case. Instead, some column names retreated back to their original names and those that did, didn't carry over the values.
Why is the spatial dataframe not correctly exporting to the attribute table in the featureclass?
Code:
import arcpy
import pandas as pd
import numpy as np
import arcgis
from arcpy.sa import *
from arcgis.gis import GIS
from arcgis.mapping import WebMap
from arcgis.features import GeoAccessor, GeoSeriesAccessor

# create spatial dataframe
sdf_AllThinningPlots = pd.DataFrame.spatial.from_featureclass("studyArea_AllThinningPlots")

# drop columns
sdf_AllThinningPlots = sdf_AllThinningPlots[['OBJECTID_1', 'Project_Na', 'Project_Nu', 'Annual_Yea', 'Treatment', 'Percent_Tr', 'Mapping_Me', 'Seeded', 'Shape_Leng', 'SHAPE']]

# rename columns
sdf_AllThinningPlots.rename(columns={
    'Project_Na': 'Project',
    'Annual_Yea': 'Year',
    'Project_Nu': 'Project_Number',
    'Percent_Tr' : 'Pct_Treated',
    'Mapping_Me' : 'Mapping_Method',
    'OBJECTID_1' : 'OBEJECTID'
}, inplace=True)

# Export to featureclass
sdf_AllThinningPlots.spatial.to_featureclass(location=r"C:\dougals_fir_reg_research\data\test_AllThinnningPlots.shp")

Spatial Dataframe:

The attribute table of the featureclass after exported to a featureclass from a spatial dataframe


texts: Arcade expression to see if current date is within two dates?

texts:
I am new to building expressions in Arcade and I am having trouble writing an expression that can be used to check start and end times of different areas and evaluate if the current date is within the start and end time for each area.  Ideally the expression would check the current time and symbolize each area based on if it is within the start and end times or not.

My layer has both a start time and end time feature in Date format.  I would like to symbolize each area as either 'opened' or 'closed' based on pre set start and end times.

Here is what I have come up with so far that doesn't quite work as I'd like:

//Expression that evaluates the status of crab catch areas

//returns open if current time is within start and end times

//returns closed if current time is outside of start and end times

var started = $feature.StartTime;

var present = Now()

IIf(started < present, 'open', 'closed')

var closed = $feature.EndTime;

IIf(closed > present, 'open', 'closed')

Right now it will not correctly symbolize the layer based on given start and end times.

Does anyone know where I messed up or another way to go about making a custom expression for this task?

This is for a map created and managed on ArcGIS Online.


texts: Opening NetCDF file in QGIS gives no extra dimensions found?

texts:
I'd like to ask if anybody here had a problem with opening nc extension file? In the beginning I've tried to install NetCDF Browser and Crayfish plugins but when I'm trying to open file with them I get a statement

  No extra dimensions found

I've found similiar topic here: Opening HDF data in QGIS on Ubuntu? but my version of GQIS (2.16.1) has HDF5 support so problem lies somwhere else.

Does anybody know where?

I'm not a native english speaker so apologise for possible mistakes.


texts: UpdateCursor: value from different rows (ArcPy)

texts:
I am new to ArcPy.

I have an attribute table with the following columns and fields

I want to fill up the column Combined as follows:
If number contains 2 categories, "combined" will be A,B.

If number contains only one then this will be pasted to "combined" column.

For example in this case first 4 rows of "Combined" will be A,B
Then 5555 is the same category so only A etc..

Update cursor uses only the same row so I cannot find a way to do this.


texts: Got Pacman instead of circle when buffering point QGIS 2.14 standalone C++ application

texts:
QgsVectorLayer *circleLayer = new QgsVectorLayer("Polygon", "Circle", "memory");

QgsFeature circleFeature = QgsFeature();

circleFeature.setGeometry(QgsGeometry::fromPoint(geometryPoint)->buffer(10000,20));

QgsVectorDataProvider *circleDataProvider = circleLayer->dataProvider();
QgsFeatureList circleFeatureList = QgsFeatureList();
circleFeatureList.append(circleFeature);

if(!circleDataProvider->addFeatures(circleFeatureList)){
    qDebug()<<"Error adding feature!";
    return;
}

circleLayer->updateExtents();

QgsSingleSymbolRendererV2 *circleLayerRenderer = new QgsSingleSymbolRendererV2(QgsSymbolV2::defaultSymbol(circleLayer->geometryType()));
circleLayer->setRendererV2(circleLayerRenderer);

QgsMapLayerRegistry::instance()->addMapLayer(circleLayer, TRUE);
mainLayerSet.push_front(circleLayer);
mpMapCanvas->setLayerSet(mainLayerSet);
mpMapCanvas->refresh();
mpMapCanvas->zoomToFullExtent();

makes Pacman (cuts off last segment of circle) instead of circle in my standalone application. Is this a way to fix this?


texts: Random sampling of many millions of Points or raster cells with ArcGIS Desktop?

texts:
I need to test the mean of two very big point shapes. One point shape has around 34.000.000 points and the other around 1.000.000. I have to do a random sample in ArcGIS 10.3.1.
Unfortunately, I haven’t found a tool to do this random sample. It should be really randomly and not for example every tenth point.
I got these points by running the tool raster to point, so if there is no solution for the point shapes, perhaps there is a straightforward way to randomly sampling a raster?


texts:
I've searched through the posts here and can't seem to figure this out. I am pulling in netCDF's from here, and processing them to display on a map. So far I have it in a GeoTIFF format that I can view with QGIS3, but it is not at the right projection/warp. I'm trying to make it so that if you look down from above it lines up on a CONUS map like this.

Here's what I'm doing so far:

export GDAL_NETCDF_BOTTOMUP=NO
gdal_translate NETCDF:OR_ABI-L2-MCMIPC-M3_G16_s20182331832322_e20182331835095_c20182331835209.nc:CMI_C01 -b 1 test.tif

And that is how I get this test.tif but is at a weird angle if you look at it in QGIS3.

So from here, I am trying to make my test.tif look like this so my GeoServer instance can slice and dice and serve the world.


texts: Setting color of a feature "street" with a FeatureRenderer

texts:
I have a layer with streets and I'm trying to change the color of the streets.  I know I asked a similar question, How do I change the colors of my layer file?,  and the answer worked for points.

This is my attempt at getting it to work with a street layer.

IFeatureRenderer featureRenderer = geolayer.Renderer;
int streetFeatureID = geolayer.FeatureClass.FindField("STREET");
IFeature feature = geolayer.FeatureClass.GetFeature(streetFeatureID);
ISymbol symbol = featureRenderer.get_SymbolByFeature(feature);
ILineSymbol linesymbol = symbol as ILineSymbol;
linesymbol.Color = SymbolCreator.getColor(255, 0, 0); /*my method which creates an RgbColor*/

This seems like it should work - there are no nulls or errors - but I'm missing a way to set the symbol.

I'd like something like this:

featureRenderer.set_SymbolByFeature(feature, linesymbol);

Since I can't find that method I think I must be doing something completely wrong.  Thanks in advance.  Using ArcEngine 10 C# VS2010


texts: WorldPop Export to Drive clip in Google Earth Engine

texts:
How I can export to drive a clipped part of a WorldPop image, I'm using the following code but crush ....

// clip fun
var clipToCol = function(image){
  return image.clip(export_geometry_2);
};

var dataset = ee.ImageCollection("WorldPop/GP/100m/pop").map(clipToCol)
                .filterDate('2017');
var population = dataset.select('population');
var populationVis = {
  min: 0.0,
  max: 50.0,
  palette: ['24126c', '1fff4f', 'd4ff50'],
};
Map.setCenter(113.643, 34.769, 7);
Map.addLayer(population, populationVis, 'Population');

Export.image.toDrive({
image: dataset.select('population'),
description: 'worldPop',
scale: 100,
region: export_geometry_2
})


texts:
I want rotate symbol, but the angle is in gradians (grade, gon) with zero on east.
How to convert this angle to degrees with zero on north ?

I try (angle*0.9 + 90), but e.g. for angle with 380 gon => 432 degrees not 72°. Any solution ?


texts: Toggling multiple toolbars on/off in QGIS

texts: In QGIS 3.12, from the View menu I go to Toolbars to toggle several of them on/off.  However, when I toggle one (either way) I'm returned immediately to the previous screen.  I've tried holding down the Control key but this action is ignored. What do I need to do?


texts:
If I import from a POST DB a set of points that are referenced by grid reference, will the long/lat decimal also be added to the file? or will I have to convert the points to long/lat myself?

In short: I am importing grid references, but want to display long/lat + store the long/lat?

Any suggestions?


texts: Reduce Vector Point Frequency in QGIS 2.18.3

texts:
I have a vector point file with 50,000+ points from a model output. I have downloaded the vector field renderer plugin resulting in a flow direction for each point. What I would like to ultimately do is have an even distribution of arrows and have more or less arrows appear when zooming in or out. I have not yet found a way to set scale ranges to the vector arrows because as a result of the plugin, it is it's own style and I do not have the ability to set a scale range.

The other thing I am having trouble with is creating an evenly spaced array of points to significantly reduce the frequency, and then apply the vector arrow style from the plugin as a work around. Is there a way to create an even point array from a rule based style setting?

Here is my problem:


texts:
i got geoserver 2.1.2 and when i wanna add the jar file ( for proxy extention for geoserver)
in the desired file "geoserver 2.1.2/webapps/geoserver/WEB-INF/lib" and clicked to "Proxy Admin Page" in the GeoServer Web UI i got this error :

HTTP ERROR: 500

class org.geoserver.web.proxy.HostnameProvider overrides final method model.(Ljava/lang/Object;)Lorg/apache/wicket/model/IModel;
RequestURI=/geoserver/web/;jsessionid=p5z1wj01ackj

Caused by:

java.lang.VerifyError: class org.geoserver.web.proxy.HostnameProvider overrides final method model.(Ljava/lang/Object;)Lorg/apache/wicket/model/IModel;
    at java.lang.ClassLoader.defineClass1(Native Method)
    at java.lang.ClassLoader.defineClassCond(Unknown Source)
    at java.lang.ClassLoader.defineClass(Unknown Source)
    at java.security.SecureClassLoader.defineClass(Unknown Source)
    at java.net.URLClassLoader.defineClass(Unknown Source)
    at java.net.URLClassLoader.access$000(Unknown Source)
    at java.net.URLClassLoader$1.run(Unknown Source)
    at java.security.AccessController.doPrivileged(Native Method)
    at java.net.URLClassLoader.findClass(Unknown Source)
    at org.mortbay.jetty.webapp.WebAppClassLoader.loadClass(WebAppClassLoader.java:366)
    at org.mortbay.jetty.webapp.WebAppClassLoader.loadClass(WebAppClassLoader.java:337)
    at java.lang.Class.getDeclaredConstructors0(Native Method)
    .

Im looking forward this error already one week and i cannot find any suggestions or solutions to handle this error :s

Need HELP please !!


texts: Merge two polygons, attributes are added - Tool?

texts:
I have a bunch of polygons (counties) that I need to merge together.  They all have various population attributes.  I need to be able to selectively merge two counties together and get a single polygon with the combined population totals for the attributes.  Is there a way to do this in ArcGIS?

I've tried the Merge tool but while that works for the boundaries, it doesn't work for the adding the two population totals together, you have to pick one or the other.  The Dissolve tool seems to operate on all polygons with a single attribute, but there is no unifying attribute that I have that determines which polygons are merged.

Thanks.


texts: ArcMap HTML Popup XSL: Images with relative path

texts:
I added a HTML Popup to a map that displays several fields data and one image that's stored in a folder next to the map.mxd file. One of the fields in the table is the filename.jpg for each image. All local.

I have a restriction, I can't use absolute paths. It needs to be relative, like "Images\P01.jpg". The problem is, if I use relative path like that, it gets related to the temporary folder where the XML file is generated and not to my MXD folder:

file:///C:/Users/demian/AppData/Local/Temp/arcCFE7/Images/P01.jpg

I see that the xsl templates have tags like <xsl:value-of select="FieldsDoc/Title" /> that gets properties from my map, but I can't find any documentation about these... Already opened lots of xsl within arcmap templates, but I couldn't figure it out...

How can I make it relative?

Perhaps there's a different approach to solve it without XSL.

I tried simply putting a html tag in the field, none worked:

For , I get
"file:///C:/Images/".

For , I get "file:///C:/Users/demian/AppData/Local/Temp/arcCFE7/Images/".

For , I get file:///C:/Users/demian/AppData/Local/Temp/Images/

Code:

<xsl:variable name="docPath" select="'C:\Users\demian\Documents\Mapa\Images\'"/>

<xsl:when test="FieldName[starts-with(., 'Imagem')]">

    <xsl:attribute name="src">
        <xsl:value-of select="$docPath"/><xsl:value-of select="FieldValue"/>
    </xsl:attribute>

</xsl:when>


texts: ArcMap MODIS NDVI for district in India?

texts:
I am working on 159 HDF files for the India subcontinent. I am trying to get min, max, standard deviation from the different districts of India using shapefiles.

How do I make all the HDF files to represent the min, max and stddev for each district using ArcMap 10.3.1?


texts: Calculating raster cell for total length of street segment with weighted-overlay on centrality value using R

texts:
I have completed the centrality calculation from the OSM database. I want to convert the data into a raster format for crossing analysis with other variables in R.
I have followed the Q&A at Convert line shapefile to raster, value=total length of lines within cell
However, I wish to add the weight value of centrality in the raster cell.
I wonder which part I could add more procedure on weight overlay calculation?
library(raster)
library(rgdal)
library(rgeos)

roads <- shapefile("edges_bom.shp")
roads <- spTransform(roads, CRS("+proj=utm +zone=37 +south +datum=WGS84"))
rs <- raster(extent(roads), crs=projection(roads))
rs[] <- 1:ncell(rs)

rsp <- rasterToPolygons(rs)
rp <- intersect(roads, rsp)
rp$length <- gLength(rp, byid=TRUE) / 1000
x <- tapply(rp$length, rp$layer, sum)
r <- raster(rs)
r[as.integer(names(x))] <- x
plot(r)

require(RColorBrewer)
spplot(r, scales = list(draw=TRUE), xlab="x", ylab="y",
       col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")),
       sp.layout=list("sp.lines", rp),
       par.settings=list(fontsize=list(text=15)), at=seq(0, 1800, 200))

library(sp)
library(raster)
library(RColorBrewer)

require(RColorBrewer)
spplot(r, scales = list(draw=TRUE), xlab="x", ylab="y",
       col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")),
       sp.layout=list("sp.lines", rp),
       par.settings=list(fontsize=list(text=15)), at=seq(0, 1800, 200))

writeRaster(r,'r.tif')


texts: Can i get earth quake data from api services?

texts:
I am working with environmental application and i am willing to add information about Earth Quake.

If Other apis like openweather,accuweather etc have services for earthquake, please let me know.

(Lets just forget about how to use or where to use.)

There are bunch of apis but they are not free.

Is there any free api services for earth quake updates?


texts: Google Maps/Earth becoming GIS solutions yet?

texts:
I've been researching internet lately to make a comparaison between Google Geo solutions ( Google Maps, Google Earth) and other GIS solutions (Free of charge or commercial, QGIS, ArcGIS, etc.). I can say I have now an idea about how they all compare, still, I wish to confirm my conclusions by whatever your experience have served you, especially of those who dealt with Google Geo in their projects.

My question is : is Google Geo (let me underline Google Earth here) becoming a GIS solution ? How are geoprocessing tools and workflows in the latest versions of Google Earth ? What does it offer in term of geo-analytical capabilities ?

Edit : I forgot to mention that I am well aware of this thread, I just need extra infos about GIS-like capabilities of Google Geo solutions.


texts: Finding Point Counts in Multiple Ring Buffers for Multiple Point IDs Using Python

texts:
I have a shapefileA that has 13 points.  Each of these points is coded with a unique ID from A to M.

I created 10 different buffers (representing different distances) around each point via the multiple ring buffer (analysis).  I now have a shapefileB that shows all of the buffers around each of the 13 points.

I have another shapefileC that contains thousands of points, where each point is associated with a unique code (A-M) that corresponds to a single point in shapefileA.  For example, a query of “ID = A” for the shapefiles will show the one point in shapefileA, the particular set of buffers in shapefileB, and the multiple points associated (that fall within the buffers) with it in shapefileC.

Here is where I need help.  What is the most efficient way, via python, to find how many points are within each buffer band for each unique ID, and then updating a new field (in shapefileB) with the number?  I know I need a cursor to go through each row, then a select by location, then a count of points, and finally a new field calculation to show the sum of points within the specific distance buffer for the specific ID.

I’m not good with nesting cursors, and it seems like this is what I’ll need.


texts: QGIS error message while running PostgreSQL code. WARNING Invalid PostgreSQL layer

texts:
I have run a code in pgAdmin it returns 43 rows of data. I took the same code ran it in QGIS and found the same row of data. However while trying to add the new data into the map QGIS comes with the error

WARNING    Invalid PostgreSQL layer

I ran two forms of this code and both came back with the same error.
SELECT pin, ST_MakeValid(parcels.geom), ST_MakeValid(address_points.geom)
FROM   parcels
JOIN   address_points
       ON ST_Contains(parcels.geom,address_points.geom)
WHERE  class = 'CLUB';

or
SELECT pin
FROM   parcels
JOIN   address_points
       ON ST_Contains(ST_MakeValid(parcels.geom), ST_MakeValid(address_points.geom))
WHERE  class = 'CLUB';

What I am trying to do is select all the parcel that contains club within, and will add the new data to a new  table (create table AS ....).
I was making the map first to see if the data was good to map, which lead to the above error.


texts: How to add delimitedtext Layer via Python in QGIS 3?

texts:
I have a problem to add a *.tab file as delimitedtext layer in QGIS 3.

uri = ("file:///{}?"
               "&delimiter={}"
               "&xField={}"
               "&yField={}"
               "&crs={}"
               "&decimal").format(tab_text,
                                        "\t",
                                        "Rechtswert",
                                        "Hochwert",
                                        "EPSG:31468",",")
vlayer = QgsVectorLayer(uri, "layer", "delimitedtext")

If I test the layer with .isValid I get a False. Did I have to name more variables?


texts:
How can you show/hide a layer (not a feature) in openlayers 5?

There is a visible parameter to the TileLayer object i'm working with, but changing that value does nothing.  https://openlayers.org/en/latest/apidoc/module-ol_layer_Tile-TileLayer.html


texts:
I am having a problem with QGIS 3.18.
I did an unsupervised Classification Using OTB tools.
Now I want to convert the result in Polygon.
But the process go up to 97 % and QGIS reports a bug.
I thought is was the size of the area, so I clipped a small size but same thing. The process goes up to 97% and Stop.
How can I get through this issue?


texts: How do I create a network from connected links

texts:
I have a bundle of links in 'some' networks. I need to find out which ones are connected to each other i.e. I need to find out which links are not connected to the 'main' network.

I have been working on a flood fill recursive method, which effectively goes from link to link, finding out which ones are connected, and recursively doing the same.  But the numbers are quite high, and my machine just bailed, so I wondered if there were other, better, or more efficient, ways of doing this.

Using arcpy

Thanks in advance


texts: Rectifying zonal statistics processing error in QGIS?

texts:
I am attempting to use the zonal statistics tool of QGIS v. 2.18.10.

I am overlaying a vector layer containing zones on a raster layer that
contains temperature data. When I use a climate raster layer from
WorlClim (http://worldclim.org/version2) everything works as expected.

However, when I use a climate raster layer from:

https://utexas.app.box.com/v/Xavier-etal-IJOC-DATA

I receive the process does not run and I receive the error message:

unsupported operand type(s) for *: 'float' and 'NoneType' See log for
more details

Where is the 'log'?


texts: ArcGIS ArcPy update attachments table rows

texts:
Some background:

We were using Collector on Android to collect data for a layer, and we found out that the pictures that were attached contained no EXIF tags at all, including the orientation tag - meaning some photos show upside down, left side up etc. Now we have to rotate every photo manually to make them appear right side up.

Apparent solution:

In order to do that I thought the following steps will be sufficient:

Extract blob data from the attachments table into .jpg files.
Edit the .jpg files and rotate them.
Extract the data from the .jpg files and update the blob field in the attachments table according to the attachment name field.

I managed to figure out steps 1 and 2, but I wonder how to get step 3 to work?

How do I update specific attachment blob fields, from .jpg files content, according to the attachment name?

Is it possible to do with a cursor?


texts: Expression for “SelectLayerByAttribute_management” with value of a cursor in the field name?

texts:
Within a script I create a new file, which has a dynamic value in the attribute name. This also changes the name of the attribute [FID] in each file. During a query in a file I have to use this attribute. So I have to create the field name for the query expression. But this is where I fail because of the syntax of all the quotation marks and backslashes.

The "normal" (undynamic) query would look like this:

arcpy.SelectLayerByAttribute_management(
    "Union",
    'NEW_SELECTION',
    "\"FID_4113_T\" <> -1 AND \"FID_4113_L\" <> -1"
)

Doc: https://desktop.arcgis.com/en/arcmap/10.3/tools/data-management-toolbox/select-layer-by-attribute.htm

My problem is that I cannot integrate the value of the cursor (row[0].zfill(4)) into the field name [FID_*row[0].zfill(4)*_T]


texts: Does QGIS Cloud support offsets for symbology?

texts:
I used the offset settings for symbols in QGIS 2.18 but when uploaded to Cloud, the offsets don't seem to exist.

Can this setting be used in QGIS Cloud or is there an alternative?


texts:
I like to think I'm relatively well-versed in datums, projections, and coordinate reference systems.

OpenStreetMap data is stored in WGS84 Lat/Lon (EPSG:4326). This CRS is geographic and therefore intended to store locations on a 3D globe. When I see this map what am I looking at?

Looking at Reykjavik it appears that the longitude lines are tending towards the North Pole at the Greenwhich Meridian (or maybe they planned the city in perspective view??) so I assume 0,90 is top-centre of the map.

Does this class as a projection, and if so which one?

Is there one projection that is commonly used for lat/lon data?

I have data stored in Oracle using EPSG:4326 and when rendered by GeoServer without a specified projection it also displays these characteristics.


texts:
I want to reproduce the dynamic US map showing the obesity percentage across the years 1984-2009 as shown in the CDC website http://www.cdc.gov/obesity/data/trends.html

I am looking for tools and code to get started on this.

EDIT 1 (Apr 22 2011)

I was able to gather the data (1994-2009) and I have uploaded it to Google docs in csv format.

Using R and googleVis I was able to generate a map for a single year (say 1994) using gvisGeoMap. But I want to create a map which will motion through years 1994 to 2009. I am not sure how to add the motion ability to gvisGeoMap. Do I need to look at some other tool?


texts: Python - Issue to create choropleth map based on GeoJSON file: AttributeError

texts:
I'm trying to build a choropleth map, I cannot figure out why I get this error:

  AttributeError: 'NoneType' object has no attribute 'get'

Here is my code:

station_nb is a dataframe with values grouped by c_ar

geop = r'arrondissements.geojson'
map.choropleth(
    geo_data =geop,
    name='choropleth',
    data = station_nb,
    columns = ['c_ar','count'],
    key_on='features.properties.c_ar',
    fill_color='YlOrRd',
    fill_opacity=0.7,
    line_opacity=0.2
)

When I check the contents of GeoJSON file key looks ok and I get the 20 expected values:

json.load(open(geop))['features'][7]['properties']['c_ar']

geop is as follows:

{'type': 'FeatureCollection',
 'features': [{'type': 'Feature',
   'geometry': {'type': 'Polygon',
    'coordinates': [[[2.351518483670821, 48.8644258050741],
.........
  'l_aroff': 'Bourse',
    'c_arinsee': 75102,
    'c_ar': 2}},

Could you help me to understand this issue?


texts: Insert query result to table (PostgreSQL)

texts:
This Query:

ALTER TABLE planet_osm_line ADD vertices FLOAT;
INSERT INTO planet_osm_line(vertices)
SELECT st_NPoints(way) FROM planet_osm_line;

does result in a empty column (NULL). What am I doing wrong?


texts: Converting Lat/Long coordinates to ESRI format

texts:
Currently I have a CSV file which has Lat/Long coordinates which I got from the batch geocoding site, http://www.findlatitudeandlongitude.com/batch-geocode/#.VkWv97fNy70.

I believe these are obtained from Google maps. I would like to map these coordinates. However, when I try to use Layer>Add Delimited Text layer, and select my Long and Lat coordinates all I get is errors.

For reference the coordinates in my CSV file are as such.

  51.494738 -9.720526

While when I upload my map into QGIS and zoom in to a region the coordinates which appear on the bottom of the screen are as such;

  152930,309827

How do I go about either converting my coordinates in my CSV so they can be read in QGIS or can I change the coordinates of my Shapefile?


texts:
I have 2 tables "dijkstra_results_geo" and "dijkstra_results_geo_aggregate".
I have a set of Linestrings in the column dijkstra_results_geo.geom_path_intermediaire (defined as geometry(Linestring, 4326)).

I would like to collect all these Linestrings into one big Multilinestring in the column dijkstra_results_geo_aggregate.path (1 row with gid=1).
I tried :

 UPDATE dijkstra_results_geo_aggregate
 SET path = (SELECT St_Union(dijkstra_results_geo.geom_path_intermediaire)
FROM dijkstra_results_geo)
 WHERE gid=1;

The query works but the "path" column (defined as geometry(Multilinestring, 4326)) remains empty.

I also tried to use St_Collect but it's empty as well.

I finally tried the query with the "path" and "geom_path_intermediaire" columns defined as "geometry", but it's the same.
Could you please help me?


texts: Polygon found with search cursor, but not update cursor

texts:
Here's my issue. I am doing two things with arcpy:

1) I search a feature class using arcpy.da.SearchCursor() and pull out a UID field and the geometry in WKB format for all polygons in the file, storing them in a dictionary.

2) I open a different, but similar feature class and use arcpy.da.UpdateCursor() to update the geometries of polygons that have the same UID as above.

The problem is not all of the polygons are getting updated, even though the exact same UIDS exist in both feature classes. Some of them update just fine, but just under half do not.

import arcpy, time

updated_rid_dict = {}
updated_rows = 0

# Grab UIDs and geometries of updated regions
with arcpy.da.SearchCursor(good_geometry_shape, field_names=['RegionID', 'SHAPE@WKB']) as cursor:
    for row in cursor:
        updated_rid_dict[row[0]] = row[1]
print "{} updated regions found".format(len(updated_rid_dict)) # returns 594

# Open an edit session and attempt to find all the regions in updated_rid_dict
edit = arcpy.da.Editor([path to my data])
edit.startEditing(False, False)
edit.startOperation()

with arcpy.da.UpdateCursor(bad_geometry_shape, field_names=['RegionID', 'SHAPE@WKB']) as cursor:
    for row in cursor:
        # See if the feature is in the updated_rid_dict
        if row[0] in updated_rid_dict:
            updated_rows += 1
            del updated_rid_dict[row[0]]

edit.stopOperation()
edit.stopEditing(True)

print "{} rows updated.".format(updated_rows) # Returns only 313

The code above doesn't even attempt to do an update. I can literally change "UpdateCursor" to "SearchCursor" and suddenly arcpy can find all the polygons in updated_rid_dict. I change it back to "UpdateCursor" and arcpy can only find some of the regions. I'm clueless as to what could be going wrong. Any ideas?

UPDATE:
Using normal cursors does not fix the problem. Further testing shows that arcpy.da.UpdateCursor() iterates through every record up to OID 16781, then stops. Does UpdateCursor have some sort of limitations I'm not aware of? 16,781 records doesn't seem like that many to me.


texts: View data of a particular year in Google Earth

texts:
For my region in India, Google Earth tiles belong to different years, 2003 - 2010. I don't mean GE's historical view feature but the imagery updates for different areas my region.

I classify Landsat data for particular types of vegetation and check out the effectiveness by overlaying the vectorized layer on Google Earth. Since I use 2010's Landsat data I want to see my layers overlaid on Google Earth's 2010's tiles when available and 2003's tiles not to appear at all. Will it be possible to do this? If yes, how can I do that?


texts:
I am using QGIS 2.1.0, and  trying to import my vector layer(parcels) with OpenStreetMap layer.
My vector layer (parcels) use EPSG 31277 coordinate system.

This is the definition for EPSG 31277:

(+proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs)

OpenStreetMap layer use  WGS 84/Pseudo Mercator EPSG 3857.

My project use WGS 84 / Pseudo Mercator EPSG 3857.
In project properties I enable 'on the flay' crs transformation.

The problem is that my parcels do not match with OpenStreetMap layer!
The difference is about 300 m. Here is a picture that shows that.


texts: Problems with adjacency matrix with poly2nb in R

texts:
I am working with a shapefile of Ohio precincts, with the hope of performing areal regression models (probably BYM2). To build the adjacency matrix, I used poly2nb from the spdep package.  As a double-check, I went back and plotted the neighbors for each county. I'm glad I did. In some counties it worked brilliantly, but in others there were numerous missing connections between precincts.
A simple example is illustrated below for Ottawa County. In the northwestern corner, you can see a precinct that lacks any connections, but which obviously borders at least five precincts.  In addition there is a cluster in the south-central area that connects among itself, but not to anything else. To say nothing of the islands.
I don't really know where to start with this, as I have no idea what would cause this.  I can manually add the adjacencies for Ottawa County, but for a larger county like Mahoning or Lucas it would be almost impossible. This is the type of thing for which there might be an obvious list of things to check for debugging/cleaning, but I don't know where to begin.


texts: Can't license ArcGIS Pro from an evaluation

texts: I have the evaluation version of ArcMap 10.3 and it is working well. However, I was trying to run the associated ArcGIS Pro 1.4 but it won't allow me in. It keeps asking me for my ESRI account details and whenever I put it, it pops out an error notification that I cannot use it as an individual. Please, I need assistance on this.


texts: Not loading all features in attribute table viewer/editor when QGIS opens?

texts:
Is there a way to have the QGIS attribute table viewer/editor not load the entire dataset when opening? I'm working with a shapefile of 100,000 grid references and this is rendering the application unusable.

I have 12GB RAM. Experiencing this behaviour on both Windows 7 and OSX.

Running latest QGIS stable version 2.18


texts: Converting lats and longs to centerpoint in ArcGIS API for JavaScript?

texts:
I'm working on a angular 2 website using the Arcgis API 3.24.  I'm fairly new to ESRI and I was wondering how I could convert a extent that has values of

spatialReference :
    {wkid: 4326}
    xmax:
        "-106.477219"
    xmin :
        "-106.728864"
    ymax:
        "35.190922"
    ymin:
        "35.014214"

to a centerpoint x and y that I could center my map on?

My code:

setSelectedItem(selItem: any) {
  const self = this;
  loadModules(['esri/geometry/Extent', 'esri/SpatialReference', 'esri/symbols/SimpleLineSymbol', 'esri/graphic']).then(([Extent, SpatialReference, SimpleLineSymbol, Graphic]) => {
    const extent = new Extent(selItem.extent[0], selItem.extent[1], selItem.extent[2], selItem.extent[3], new SpatialReference({wkid: 4326}));
    self.isDrawToolVisible = selItem.useExtent;
    self.selectedItem = selItem;
    const color: any = '#004065';
    const symbol = new SimpleLineSymbol(SimpleLineSymbol.STYLE_SOLID, color, 4);
    self.extentLayer.clear();
    self.extentLayer.add(new Graphic(extent, symbol));
    if (self.cartMap.graphics) {
      self.cartMap.graphics.clear();
    }
    if (selItem.useExtent && selItem.clipGraphic) {
      self.cartMap.graphics.Add(selItem.clipGraphic);
    }

    self.downloadFinishedItems = self.cartService.getItemsFinishedDownload();
  });
}

UPDATE :

When I try to center the x and y return with NaN:


texts: Rotating only one feature of a point layer

texts: I am working on a gas network utility and I create some point layers in order to insert some symbols of special forms. However, each symbol should be attached the network with a specific angle in order to fit it with the network. Rotating actions applies to all point symbols of the layer. Is it possible to rotate only one feature?


texts: Modify the Zoom Buttons in ArcGIS.com JavaScript Application

texts:
I'm writing a JavaScript application to view web maps from ArcGIS.com.  It is similar to the sample found on Esri's Samples page:

http://developers.arcgis.com/javascript/samples/ags_createwebmapid/

I've run the Chrome Developer Tools on that page and it looks like the Zoom buttons in the upper left are autogenerated by the referenced JS libraries.

Right now I just want to add a "Full Extent" button, but I may make additional changes later on.  Is there a straightforward way to modify those buttons in the JavaScript code?


texts:
I have a one question: it's possibile to create an offline map of openrailwaymap?

I would create an offline map using openstreetmap and after add the openrailwaymap layer...


texts:
Is there an easy way to restore databases from PGadmin 3 into Pgadmin 4? Right now data on both, and I want to gather all my data on the same server. The port on PGadmin 3 is 5432 and 5433 on PGadmin 4. My PostgreSQL versions is 10.1 on pgadmin 3 and 11.3 on pgadmin 4.

Previously, with help I succeeded in doing so by using Notepad++, but I can't figure out how we did then.


texts: Placing point symbol on center of polygon in ArcMap?

texts:
I'd like to symbolize polygons with a point at the centroid. Does anyone have suggestions how to accomplish this without creating a new shapefile?

I'm using ArcGIS 10.0.


texts: Counting Characters using Field Calculator - ArcMap

texts:
I am trying to use field calculator to populate a column based on the the number of characters in another column. I am using the Python Parser. I am not sure if the .count() function is appropriate in this case or if I am even using it right.

My pre-logic script is:

def reclass(f1):
     if f1.count() >= 4:
         return "Yes"
      else:
         return "No"

Label = reclass( !HUB! )


texts: Calculating latitude and longitude values for each pixel corner?

texts:
I have a TIFF image file and I need to extract the latitude and longitude values (in decimal degrees) for each pixel's corners. I have the latitude and longitude values for the top right, top left , bottom right and bottom left corners. Along with this, I also have the pixel values for each pixel in the image.

How do I perform this task in ENVI or in IDL?

Do I need to some information about the TIFF file in order to proceed?


texts:
I'm having a hard time figuring out the MBG orientation and how to represent it in a accurate way. My work consists in mapping boulder deposits (shapefiles) and with the MBG convex hull to assign them a direction. However I don't understand how to interpret the data on the MBG_Orientation column? I've understood it's on 180 degrees only but which direction means which number (north, south, east & west- wise)? For example most of them are around 80 & 90 degrees which I would assume is a eastern direction but I'm not sure and would like to understand how it works.
Also, I would like to use the rotation symbology so as to show direction with arrows...and this gets event more complicated for me, if you could help!


texts: Convert geometry to latitude-longitude using PostGIS ST_Transform

texts: My table has a geometric field, in EPSG:3857, of type point. I need the latitude and longitude, in EPSG: 4326, of each point.  Could someone give me a valid example of how to use ST_Transform().


texts: I was wondering if it is possible to display a map in ol3, which is constructed as WCS in MapServer/GeoServer?


texts: Counting number of NULL values in each row using QGIS?

texts:
I'm new to QGIS and am trying to figure out how to count the number of NULL in each row and reflect the number in a new field.

How would the expression look like?

The image is what the table looks like.


texts: How to get a search cursor on a particular version in arcpy

texts:
I have a versioned database in ArcGIS and I am wondering how can I tell arcpy what version to query. The closest I have got to an answer has been building a particular connection file that specifies the version to connect to, or use ChangeVersion_management which only works on layers and views.

I can easily do that in a published map on the web using ArcGIS JS API by setting the version using layer.setGDBVersion(version)


texts: QGIS Segmentation fault when manipulating QgsGeometry

texts:
So I am running QGIS 2.8.1 on Ubuntu 14.04. This following problem is happening in my script, but also occurs directly in the QGIS python console:

lay = iface.activeLayer()
lay.featureCount()
1L
lay.type()
0
prov = lay.dataProvider()
prov.geometryType()
3
geom = lay.getFeatures().next().geometry()
geom
<qgis._core.QgsGeometry object at 0x7fff4d5f90e8>

Next, if I do pretty much anything with that QgsGeometry object I get a dizzying array of Segmentation faults. I noticed it first when I was trying to intersection tests with other geometries, but even simple calls like geom.type() results in :

Program received signal SIGSEGV, Segmentation fault.
0x00007fffefe0c3ca in _int_free (av=0x0, p=0x5042970, have_lock=0)
at malloc.c:3943

but mostly it has to do with the GEOS related functions:

Program received signal SIGSEGV, Segmentation fault.
0x00007ffff123798f in GEOSGeom_destroy_r () from /usr/lib/libgeos_c.so.1

is this a known bug? http://hub.qgis.org/issues/777


texts: I have a SHP file of sub-areas of U.S. counties (roughly about the size of a Census tract) that I've reprojected into a 'true distance' projection (World Equidistant Cylindrical) in QGIS 1.8. However, when I measure the length of the sub-area border segments I get results that don't make sense (3 meters for instance). I've done this before in earlier versions of QGIS and my results seemed fine. Any idea what the problem could be? Thanks for your help.


texts:
I can finish drawing a feature with a double click or F2, but I cannot run a geoprocessing tool with the same ease.

I tried looking to Esri for guidance but the page itself only makes a reference to a "Mouse Click" to run a tool, but no shortcut key...

  Run a geoprocessing tool by clicking the Run button at the bottom
  of the Geoprocessing pane.

https://pro.arcgis.com/en/pro-app/help/analysis/geoprocessing/basics/run-geoprocessing-tools.htm

Using ArcGIS Pro, how might I make a custom key, if a shortcut key for this does not already exist?


texts: how to calculate cumulative cost of movement, accounting for wind direction, in R?

texts:
I am trying to reproduce the calculation of the duration of sail-powered navigation as described in an academic article (http://www.tandfonline.com/doi/full/10.1080/15230406.2017.1403376), but using R instead of ArcGIS.

Following this nice guide (https://mgimond.github.io/Spatial/raster-operations-in-r.html), I am planning to use the R package ‘gdistance’ to calculate cumulative distances, but I am stuck in the very early steps of the procedure. Following the mentioned website, first one should have a cost raster, which in my case would be the inverse of the wind speed (i.e., the pace of the navigation speed; let’s say hours/km). The ‘gdistance’ package (https://cran.r-project.org/web/packages/gdistance/vignettes/gdistance1.pdf) allows to calculate:

(1) a transition matrix (using the ‘transition()’ command) that
defines a traversing cost,
(2) a geocorrection that adjusts the cost for 'true' local
distance, and
(3) an accumulated cost.

What is not clear to me is how an horizontal factor (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-the-horizonal-and-vertical-factors-affect-path-distance.htm) can be put to work here.

In essence, the ideal procedure should be able to take into account a wind direction raster and to increase or reduce the cost of traveling (i.e., the sailing time) according to whether or not the movement along adjacent cells is (for instance) against or in the same direction of the wind. In other words, what is not clear to me is how a second cost, different from a wind pace raster, can be taken into account using ‘gdistance’.


texts: How to overwrite exit event when using PyQGIS?

texts: I'm developing a plug-in using PyQGIS and when I click on the exit button the QGIS executes some code defined by itself (it asks if I want to save the changes and then exit). I'd like to overwrite this executing a code defined by me, but I can't find how to do it. How can I do this?


texts: Is there any way to reset feature cursor?

texts:
When we search within a feature class, it returns a feature cursor. To get features from a feature cursor, I use following code snippet:

private List<IFeature> GetFeatureList(IFeatureCursor featureCursor)
{
    var features = new List<IFeature>();
    IFeature feature;
    while ((feature = featureCursor.NextFeature()) != null)
        features.Add(feature);
    return features;
}

But problem is for some interfaces and methods, we need this featurecursor. But before passing this feature cursor to that method, I have already traverse it. So, using IFeatureCursor.NextFeature method, I don't get any other feature here. Is there any method or interface or anything that can reset the feature cursor and make the feature cursor re-usable?


texts: Converting MGA coordinates for Google Earth?

texts:
I am new to looking up coordinates, etc.

However, I have some coordinates that I want to look for, but they are in MGA.

I am wanting to look up these coordinates on Google Earth and then try and export them to my vms.

From what I can discover Google Earth does not like these MGA coordinates and I have zero know how on trying to convert them to a readable coordinate to allow me to do the above.

Is anyone able to help?


texts:
I am trying to run a code in R but I get an error. I am conducting a logistic regression in order to predict a raster by another raster. the code is as follows:

> sinks <- raster("rasterized.tif")
> tpi <- raster("TPI.tif")
> model <- glm(sinks[]~tpi[], family= binomial)
> summary(model)

Call:
glm(formula = sinks[] ~ tpi[], family = binomial)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-0.5368  -0.3780  -0.3774  -0.3766   2.3877

Coefficients:
             Estimate Std. Error   z value Pr(>|z|)
(Intercept) -2.606015   0.001022 -2551.117   <2e-16 ***
tpi[]       -0.321924   0.038967    -8.261   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7496150  on 14968071  degrees of freedom
Residual deviance: 7496082  on 14968070  degrees of freedom
  (15480 observations deleted due to missingness)
AIC: 7496086

Number of Fisher Scoring iterations: 5}

But then I try to use the model for another raster(tpi_test), but I get the following error:

>test <- raster("tpi_test.tif")
> predict(test, model, type="response")
Error in p[-naind, ] <- predv :
  number of items to replace is not a multiple of replacement length
In addition: Warning message:
'newdata' had 690549 rows but variables found have 14983552 rows

I am very new to R. Can anybody tell me what the problem is?


texts:
I'm playing around with a public NTRIP caster but I'm very new to all this. I'm using BKG Ntrip Client to read the stream and the file I get has around 45 lines for each sample. Each line begins with the stream code, then G01, G06, etc. Does that mean the caster can view 45 satellites? I read somewhere that usually only at most 10-15 satellites are visible.

EDIT: Here's the first 2 lines I got:

ROBI4 G01 C1C   23698316.588 L1C  124535486.180   15 S1C   40.000 C2W   23698320.484 L2W   97040899.394   15 S2W   22.000 C2X   23698320.966 L2X   97040571.418   15 S2X   39.000 C5X   23698321.324 L5X   92997563.294   15 S5X   46.000
ROBI4 G06 C1C   25005922.799 L1C  131407234.494   13 S1C   34.000 C2W   25005927.678 L2W  102395247.228   13 S2W   17.000 C2X   25005928.446 L2X  102394614.238   15 S2X   35.000 C5X   25005928.232 L5X   98128172.506   15 S5X   40.000


texts: Performing Point Distance analysis using Basic level license of ArcGIS Desktop?

texts:
I am using ArcView (ArcGIS Desktop Basic) 10.1 and I need to perform a point distance analysis.

What are the steps to determine the distance from A to B-Z, B to A-Z, C to A-Z, etc?


texts: Mapping shapefile attributes to raw data

texts:
I have raw address data stored in a sql server database (2012) which comprises address, city/suburb, state and generally country (Australia).  I would like to take this address data and following geocoding to long/lat coordinates, match these points to a respective shapefile containing administrative boundary data/polygons.  My aim is to add columns to my address data from attributes from the shapefile e.g. LGA_NAME for the purpose of mapping within reporting tools like report builder 3 or Tableau.  I have a copy of QGIS I am experimenting with but haven't found a way to do what I want as yet.

I suspect this would be a fairly straightforward task and the point in polygon concept suggests it should be possible.

I have looked for a suitable tutorial for this but am yet to find anything.

Can anyone help?


texts:
"Strong OS Stack Knowledge" is listed as a required skill on a job specification I am interested in. I haven't come across this term before. It may be important to note that I live in the UK and I think OS may either mean Operating System, Open Source or Ordnance Survey.

Any ideas?


texts:
Background
I'm trying to create parallel lines that go through each parallel points. In this points arrangement, it's North-South lines:

I use Python and the libraries: pygeos, geopandas, and shapely
Difficulties
My current approach is to create Delaunay lines, then filter it based on its azimuth. I keep the deviation around +- 10 degree from 0 and 180 (0 & 180 based on minimum_rotated_rectagle orientation)
Currently here's my result:

Here I see there're two things that should be done.

Connect the lines
Connect points that are not in any lines

To solve both of this problem, I'm thinking to create a connecting line for each line/point to nearby lines (filtered by threshold distance (x and y); ex: line 8  has nearby lines -> 6, 7, 9). When connecting the lines the connecting line created must be either at the start or end of the line.
Using the created connecting lines, then I filter it again using the azimuth.
Questions

What line should I use, is it hausdorff lines? (if it's, is there library that create the line? not only calculate the distance).
Is there any better approach for this problem?

I've tried to use voronoi and buffer to simulate neighboring points, however it's not possible because there could be deviation in points position, therefore not creating a straight line.
Desired Output
Here's the desired output. The lines connected to points in its parallel.

For point in green circle, it will be combined to neighboring line based on the shortest distance.


texts: Which data format can be used with both QGIS Desktop and Android?

texts:
My data is in MS Access geodatabase(creating with ArcGIS). I want use Qgis desktop and Qgis android (and dropbox)

Is it required to convert my mdb file to another format? (Qgis on Android can't read my mdb file)
Which format should I use?


texts: What is the difference between 10 metre U component of wind and U component of wind

texts:
I am trying to read GRIB files which have a number of messages including

10 metre u component of wind
10 metre v component of wind
u component of wind
v component of wind

I cannot seem to find any mention of a u component of wind on the internet, only 10 metre.

What is the difference?


texts: I would like to use the stepwise regression method to compute time series analysis in GEE, is there a way to implement this regression model in GEE.


texts: Draw Triangle Polygon with Dynamic Display

texts:
I'm trying to draw a custom polygon in Dynamic Display but its really not working at all. I can get it to work if I do this:

IPointCollection pc = seg.setRectangle( mapcontrol.map.extent.envelope) as IPointCollection;

But I want a custom polygon.  I tried this following code but I'm getting -2147220968

       IPoint p0 = new PointClass();
        IPoint p1 = new PointClass();
        IPoint p2 = new PointClass();
        p0.SpatialReference = gcs;
        p1.SpatialReference = gcs;
        p2.SpatialReference = gcs;

        p0.PutCoords(-94.569168d, 46.697493d);
        p1.PutCoords(-94.052811d, 46.711619d);
        p2.PutCoords(-94.51561d, 46.520131d);

        p0.Project(mapSP);
        p1.Project(mapSP);
        p2.Project(mapSP);

        ILine line = new LineClass();
        ILine line2 = new LineClass();
        ILine line3 = new LineClass();

        line.SpatialReference = mapSP;
        line2.SpatialReference = mapSP;
        line3.SpatialReference = mapSP;

        line.FromPoint= p0;
        line.ToPoint = p1;

        line2.FromPoint= p1;
        line2.ToPoint = p2;

        line3.FromPoint= p2;
        line3.ToPoint = p0;

        ISegmentCollection seg = new PolygonClass();
        seg.AddSegment(line as ISegment);
        seg.AddSegment(line2 as ISegment);
        seg.AddSegment(line3 as ISegment);
        pc = seg as IPointCollection;

                           dynamicSymbolProps.set_DynamicGlyph(esriDynamicSymbolType.esriDSymbolFill, fillGlyph);
try
{
  dynamicDisplay.DrawPolygon(pc);
}
catch (Exception ex)
{
  string m= ex.Message;
}


texts:
I am attempting to remove duplicates in a feature class that are generated by a Python script.  At the end of my script I have the code block below; I haven't implemented it yet, but would like to see if this is the correct way to use arcpy.da.Editor as the errors that I have been receiving are in regards to being outside of an edit session.

Is arcpy.da.Editor being used correctly here?

`

        fields = ['OBJECTID', 'SRNumber', 'ItemLoc_9']

        table_rows = []
        delete_rows = []

        where_clause =  """ ItemLoc_9 = ' '"""

        with arcpy.da.Editor(appendClass) as edit:
            with arcpy.da.SearchCursor(appendClass, fields, where_clause) as cursor:
                for row in cursor:
                    if [row[1], row[2]] in table_rows:
                        delete_rows.append(row[0])
                    else:
                        table_rows.append([row[1], row[2]])

                del table_rows

        with arcpy.da.UpdateCursor(appendClass, ['OBJECTID']) as cursor2:
            for row in cursor2:
                if row[0] in delete_rows:
                    print 'Deleting record: OBJECTID = ' + str(row[0])
                    cursor2.deleteRow()

                    print '\ndone'`


texts:
I'd like some general advice on creating dynamic maps. I have a bunch of geospatial points saved on SQL Server and PostgreSQL. People will enter new data into the databases from time to time. I want those points to be automatically plotted into a series of base layers I have already created with QGIS. The plotting does not have to be instantaneous, but it would be nice if the map was updated every day. I've got a few leaflet maps that I especially like, but I'm not tied to them. I have a decent bit of SQL, some python, marginal amounts of HTML and PHP experience, and practically no JavaScript. To be perfectly honest, I don't know where to start in figuring out how to this. I have very little as far as budget is concerned, but lots of free time to learn how to program better. I have been doing some research and I'm still pretty confused on where to start. I was wondering if anyone had a good resource on where to get started.

Right now, I am thinking the best thing for me to do would be learning geoalchemy2 on python to create the type of map I want to create. I'm pretty sure geoalchemy2 can do this, but could anyone confirm?


texts: Is it possible to sync multiple radio buttons on a QGIS 2.0 custom feature form?

texts: Is it possible to sync multiple radio buttons on a QGIS 2.0 custom feature form to one attribute? For example i have three choices (owner, renter, and other) and choosing the radio button would fill the same field with the button of choice.


texts: Is there a 'find' function in Openlayers?

texts:
Is there a method in Openlayers, where you can search a layer's attributes and zoom to the result, similar to GIS or Google Maps?

Eg. if you had a layer of streets, you could type in a street name, which would search through a 'streetname' field and zoom/highlight the result in your web map.


texts: Making grayscale map from multiple band satellite image in ArcGIS for Desktop?

texts:
I have satellite multiple band map, I needs to make grayscale map,what should i do?
what is the steps?
which band and format of input are to be used in this?
How to make attribute table of this?

In my subject i want to make classification or extraction of some pixel value of grayscale map.for that i requires attribute table. please say me about this all steps. I am new learner in this.


texts:
I am studying how the accesibility and mobility is affected by multiple factors such as traffic, extreme water events, etc. I am looking to add cost barriers, some of them should affect the speed whereas some other should block some road parts.

My goal is to show the difference between different scenarios.

The data I have includes the streets and their speed and some information for traffic and flooding areas (bluespots).

How can I add different type of barriers?


texts:
I need to normalize the values attained from a band ratio to 0-255.
I'm using QGIS. How do I do this?
Any ideas?


texts: Understanding resource.xml of ArcGIS Network Analyst?

texts:
I want to translate the direction result of network analyst to turkish. I surf in the interent and found i have to change resource.xml file . the file located in arcgis install directory

...\ArcGIS\Desktop10.0\NetworkAnalyst\Directions

I want to find a help file or dictionary file to translate each tag and attribute of the resource.xml file.

For example what does monobranch means in the xml file?

Is there any help file to understand this file?


texts: When creating annotation from labels is it possible to feature-link another annotation layer?

texts:
I have an annotation layer that is currently being used in the flex viewer and I would like to add a secondary annotation layer/line so that it sits directly underneath the first.

I am able to link the new annotation layer to the polygon I am creating it from, but I wonder if there is a way to also link it to the annotation layer I already have; as I am having trouble getting the annotation to sit underneath the first one.

I tried setting the offset, but as it offsets from either the center or the edge of the block, the new annotation does not line up within the top right corner underneath the original annotation. Editing is off the table as I would have some 30,000 records to manipulate manually.

Is there any way to link these labels to the original annotation layer before I turn them into annotation?


texts:
So this seems to be an issue because I've seen several questions concerning errors when merging or mosaicking rasters in r, including mine. Sometimes I have no issues merging rasters and others I get the follwing.
Error in compareRaster(x, extent = FALSE, rowcol = FALSE, orig = TRUE,  : different origin

I'm confused because I can go into ArcMap and do the same thing with no issues or I can figure out the differences in tolerance Raster merge attempt, different origin error and set the tolerance accordingly to overcome the error.
So if ArcMap ignores it completely and I can just overcome the error by changing the tolerance what is the purpose of the tolerance at all? Are there any negative effects from increasing the tolerance or can I just put it at the maximum and ignore it altogether?
I couldn't find anything referring directly to what the purpose of the tolerance is and why it is important.


texts: I have a point file saved as a shapefile, now I want to produce a heatmap of the values. However the heatmap is projected upside down. I suspect my y coordinates are inverted or negative but I do not know how to fix it. Can somebody help?


texts: Is there a way to run a Python script inside pgAdmin 4 without relying on e.g. the QGIS Python console to do so?

texts:
I have a python script that contains within it an SQL command. When I execute this through QGIS Python console it also runs smoothly and returns the expected results without any problems.

I would like to know whether it's possible to run a Python script in pgAdmin? Is there a possibility to add like a python console into pgAdmin or something else to achieve this?

The reason I ask is simply because if possible I would like to run everything through one programme in this case pgAdmin as the bulk of my steps are in SQL and it's just one step where I am using Python.

It just would be really helpful to know if it's possible and how to do it or if it would be more effort than it's worth since the processes as they are now all work.


texts:
I created a page where people can upload a gpx or kml file. Openlayers makes this a piece of cake.

But now I would like to provide some kind of editing to select a piece of the multilinestring.

What I was thinking

people can select a piece by putting two points on the route. "The system" copies the selected piece (between the two points) to a new source
people can trace a route. They click on the route to set a starting point, they move along the route (with snapping) and a final click to set a ending point.
people click to set a starting point, they move along the route (with snapping) and click to add a point on a turn and double click to end the draw

I have no problems with loading the file, adding snapping, drawing and editing. I'm looking for tips on a strategy to do the final part.

The last idea is probable the simpliest but i want to avoid that one (and probably the second idea for the same reason). In case of long routes it would take a lot of time and i don't want people giving up because it takes to long or is to much work.

I was wondering what other ideas are out there. Some code would be good but a push in the right direction is ok.


texts:
I want to access a couple of GPX files online to extract the tracks and waypoints for some feature classes.  When I download the files, ogrinfo (version 2.4.1 on a Mac)works fine.

ogrinfo /Users/ianhorn/Downloads/KAT_Motorcycle.gpx
INFO: Open of `/Users/ianhorn/Downloads/KAT_Motorcycle.gpx'
      using driver `GPX' successful.
1: waypoints (Point)
2: routes (Line String)
3: tracks (Multi Line String)
4: route_points (Point)
5: track_points (Point)

However, I would prefer to just grab the URL and process that way using vsicurl, but ogr spits out an error saying it cannot open the datasource with the driver.

ogrinfo /vsicurl/http://annessky.net/kat/gpx/KAT_Motorcycle.gpx
FAILURE:
Unable to open datasource `/vsicurl/http://annessky.net/kat/gpx/KAT_Motorcycle.gpx' with the following drivers.

Is it not possible to do this?

Update.  Solution provided by @Maptastik using python, https://repl.it/@maptastik/katmotorcyclegpxdownload


texts: Generate polygons of city blocks from streets with higher accuracy

texts: I have a model on QGIS that enumerate all the blocks within 1km radius with the use of Polygonize tool and OSM street data (downloaded from Geofabrik). However, the space between some wide streets/avenues are being counted as blocks. I painted some of them blue for reference. I tried to use Dissolve tool to remove these or join them with the closest 'neighbour', but it didn't work.


texts: Avoiding import arcpy error (RuntimeError: NotInitialized) when using IDLE?

texts:
ArcGIS Version: ArcGIS 10, Python 2.6.5
Machine: Windows 7 64bit (Intel)

I'm receiving the following error when importing arcpy using the the IDLE (gui and cmd prompt) and when using pythonwin.

I had removed python26 and python25 from my system and have done a fresh install of ArcGIS 10 and python. Therefore am calling IDLE from here: "C:\Python26\ArcGIS10.0\pythonw.exe"

The Error:

**IDLE 2.6.5

      import arcpy
      Traceback (most recent call last):
        File "", line 1, in
          import arcpy
        File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy__init__.py", line 17, in
          from geoprocessing import gp
        File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing__init__.py", line 14, in
          from _base import *
        File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing_base.py", line 568, in
          env = GPEnvironments(gp)
        File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing_base.py", line 565, in GPEnvironments
          return GPEnvironment(geoprocessor)
        File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing_base.py", line 521, in init
          self._refresh()
        File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing_base.py", line 523, in _refresh
          envset = (set(env for env in self._gp.listEnvironments()))
      RuntimeError: NotInitialized**

EsriUK and Esri tech support have recommended fresh reinstall and registry edits = but this has not worked. I'm waiting on them for a solution.

I've investigated whether it was a module arcpy path issue. There are a couple of links on the Esri forums and on the ArcGIS help regarding paths. Also, a possible related thread on gis.stackexchange.

I followed the advise from 'kimo' on the Esri forum thread. But still the import error.

My system variable (that I've edited) is:
PYTHONPATH = "c:\python26";"C:\Python26\ArcGIS10.0";"C:\Python26\ArcGIS10.0\Lib\site-packages";"C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy"
I have checked I have the contents of "C:\Python26\ArcGIS10.0\Lib\site-packages\Desktop10.pth". The contents are:
C:\Program Files (x86)\ArcGIS\Desktop10.0\bin
C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy
C:\Program Files (x86)\ArcGIS\Desktop10.0\ArcToolbox\Scripts

However, the traceback message shows that a .py script is found in the folder: C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy. Therefore, suggests the arcpy module is found by the interpreter but then not initialized(?).

Arcpy can be called from the Python Window within ArcMap and ArcCatalog okay. It is only when called from a Python interpreter that the 'import arcpy' fails.  Therefore, arcpy is initalized okay by ArcMap/ArcCatalog, but not IDLE - how can this be?

I previously had an install of python25, then python26 then after that upgraded to ArcGIS 10. Therefore, there was quite a lot of legacy python on my machine. However, I've uninstalled all and registry cleaned using ccleaner.

Can anyone suggest another line of enquiry as I'm not experienced with python enough to solve initialize/runtime errors?  I'd like to avoid reinstalling the OS too!


texts: Getting LonLat of a specific location in a specific projection on osm / google maps?

texts: How can I get the coordinates in, for example EPSG:4326 of a specific location, in either OpenStreetMaps or Google Maps?


texts:
I am trying to connect to my Google Maps Engine Map using WMS in my leaflet application. I am able to connect just fine with ArcMap.

However when using leaflet, I get 'no data for requested layers'. Here is a jsfiddle.

// create a map
map = new L.Map('mymap');

// set the map's starting view
map.setView( new L.LatLng(22,79), 4 );

L.tileLayer.wms("https://mapsengine.google.com/15658084116283052074-10602818151822182906-4/wms/", {
  layers: 'x',
  format: 'image/png',
  version: '1.1.0',
  transparent: true,
  attribution: "",
  tiled:true
}).addTo(map);

Do I have the wrong settings/layer?

Edit: GetCapabilities returns...

<Layer queryable="0">
    <Name>15658084116283052074-10602818151822182906-4</Name>
    <Title>x</Title>
    <Abstract/>
    <CRS>EPSG:3857</CRS>
    <CRS>EPSG:3785</CRS>
    <CRS>EPSG:900913</CRS>
    <Layer queryable="0">
        <Name>15658084116283052074-00933023778830223306-4</Name>
        <Title>x</Title>
        <Abstract>Raster layer</Abstract>
        <CRS>EPSG:4326</CRS>
        <EX_GeographicBoundingBox>
          <westBoundLongitude>-100.00045742818516</westBoundLongitude>
          <eastBoundLongitude>-92.99988641904172</eastBoundLongitude>
          <southBoundLatitude>39.99928465598993</southBoundLatitude>
          <northBoundLatitude>47.00075398041749</northBoundLatitude>
        </EX_GeographicBoundingBox>
    </Layer>
</Layer>


texts:
I use the "statistics by category" tool from QGIS to learn what grain yields I got from different combinations of scenarios like:

Soil A + Hybrid A + Planting Date A = Yield n
Soil B + Hybrid A + Planting Date A = Yield n2 ... so on and so forth.

Here`s what I do today:

I take the yield data from the combines that come as points and convert to SHP.

The shapefile now have 1 column only (Yield).

Then I use the "Join attributes by location" tool to extract data from other polygon layers that I want to analyse (like Planting maps with Hybrids locations, soil type maps etc) and add those info into my yield shapefile.

Then I run the statistics from this shapefile that has each point with all the atributes that I want to analyse.

It does work but it`s a really time consuming process, specially because my yield shapefile contains millions of points so every step takes several minutes to be processed.

Any Ideas on a more efficient procedure?

Here`s a screenshot of a final product that I have today and that I would be looking for something similar with a different solution.


texts: Can I run tools from QGIS using shell commands?

texts: Is it possible to script rutines in QGIS (python) and then runing them remotely? To illustrate, can I through a shell (ssh or telnet) ask to QGIS to transform  a shapefile into a gjson format? Mostly, I'm thinking in the interaction of QGIS with another software like R for example, where I can generate command lines. (I'm not talking about QGIS internal command shell)


texts: When i open some page, i want to show only WMS layers. When clicking on some parcel, as response, it should call WFS, and show information of that parcel, and only that one. How can i manage that with Openlayers 2.13?


texts: When I use smaller fonts in the legend (anything less than 20pts) there's an awkward gap between the items. I have set all the relevant spacing options in the Legend Properties (i.e. space between symbols, symbol label space, and line space) to 0, but the gap remains. Am I missing something?


texts: replace GIS object with illustrator symbol

texts:
My workflow entails exporting ArcGIS Layers into Adobe Illustrator for cartographic output purposes. One problem I have at the moment is that I cannot replace my points with an Illustrator symbol. For example, points from ArcGIS are imported into Adobe Illustrator as circles. I cannot replace a circle with an Illustrator-symbol. What I can do however, is replace a single anchor point with a symbol. Is there a way to convert a circle to a single anchor point? or are there other ways in which I can replace GIS objects with Symbols?

So my question is somewhat illustrator based, but I thought as this is such a common GIS workflow, is may be answered here.


texts: The underlying connection was closed: The connection was closed unexpectedly, WCF and ESRI Objects

texts:
I am getting this error when using ESRI objects like Geometry,  Envelope etc in a WCF service. I have a class that some properties (string, int, double) but when I add ESRI type as part of this class, it throws this error. Without ESRI objects, it wroks fine. I have tried to change the settings for timeout but that didn’t make any difference. I am hoping to get some feedback from other folks who have implementes ESRI objects from WCF service.

The underlying connection was closed: The connection was closed unexpectedly.

System.ServiceModel.CommunicationException was unhandled by user code
  Message=The underlying connection was closed: The connection was closed unexpectedly.
  Source=mscorlib
  StackTrace:
    Server stack trace:
       at System.ServiceModel.Channels.HttpChannelUtilities.ProcessGetResponseWebException(WebException webException, HttpWebRequest request, HttpAbortReason abortReason)
       at System.ServiceModel.Channels.HttpChannelFactory.HttpRequestChannel.HttpChannelRequest.WaitForReply(TimeSpan timeout)
       at System.ServiceModel.Channels.RequestChannel.Request(Message message, TimeSpan timeout)
       at System.ServiceModel.Dispatcher.RequestChannelBinder.Request(Message message, TimeSpan timeout)
       at System.ServiceModel.Channels.ServiceChannel.Call(String action, Boolean oneway, ProxyOperationRuntime operation, Object[] ins, Object[] outs, TimeSpan timeout)
       at System.ServiceModel.Channels.ServiceChannel.Call(String action, Boolean oneway, ProxyOperationRuntime operation, Object[] ins, Object[] outs)
       at System.ServiceModel.Channels.ServiceChannelProxy.InvokeService(IMethodCallMessage methodCall, ProxyOperationRuntime operation)
       at System.ServiceModel.Channels.ServiceChannelProxy.Invoke(IMessage message)
    Exception rethrown at [0]:
       at System.Runtime.Remoting.Proxies.RealProxy.HandleReturnMessage(IMessage reqMsg, IMessage retMsg)
       at System.Runtime.Remoting.Proxies.RealProxy.PrivateInvoke(MessageData& msgData, Int32 type)
       at TEST.Services.Test.BoundaryService.IBoundaryService.GetBoundarySetSpatial(String boundarySetName, String boundaryName)
       at TEST.Services.Test.BoundaryService.BoundaryServiceClient.GetBoundarySetSpatial(String boundarySetName, String boundaryName) in C:\Projects\TESTGeoServices\TEST.Services.Test\Service References\BoundaryService\Reference.cs:line 57
       at TEST.Services.Test.BoundaryServiceTest.BoundaryService_GetBoundarySetSpatial_Test() in C:\Projects\TESTGeoServices\TEST.Services.Test\BoundaryServiceTest.cs:line 26
  InnerException: System.Net.WebException
       Message=The underlying connection was closed: The connection was closed unexpectedly.
       Source=System
       StackTrace:
            at System.Net.HttpWebRequest.GetResponse()
            at System.ServiceModel.Channels.HttpChannelFactory.HttpRequestChannel.HttpChannelRequest.WaitForReply(TimeSpan timeout)
       InnerException:


texts: Specifying map service parameters when publishing with ArcPy?

texts:
I publish map services on a ArcGIS server using the typical arcpy script below. I need to better specify some service parameters. For example, the service should be available in WMS instead of KML, with anti-aliasing and some specific pooling parameters.

Is there a way to specify these parameters using some arcpy functions?

connection6546543 = 'C:/blabla/arcgis on myserver.eu (admin).ags'
folder1116 = "Myfolder"

def publishService(mxdPath, service_name, summary, tags):
    sddraft = 'C:/services/' + service_name + '.sddraft'
    if os.path.isfile(sddraft): os.remove(sddraft)

    # Create service definition draft
    map = arcpy.mapping.MapDocument(mxdPath)
    arcpy.mapping.CreateMapSDDraft(map, sddraft, service_name, 'ARCGIS_SERVER', connection6546543, True, folder1116, summary, tags)

    # Analyze the service definition draft
    analysis = arcpy.mapping.AnalyzeForSD(sddraft)

    # Print errors, warnings, and messages returned from the analysis
    print "The following information was returned during analysis of the MXD:"
    for key in ('messages', 'warnings', 'errors'):
        print '----' + key.upper() + '---'
        vars = analysis[key]
        for ((message, code), layerlist) in vars.iteritems():
            print '    ', message, ' (CODE %i)' % code
            print '       applies to:',
            for layer in layerlist:
                print layer.name,
            print

    # Stage and upload the service if the sddraft analysis did not contain errors
    if analysis['errors'] == {}:
        print 'Create service definition...'
        sd = 'E:/services/' + service_name + '.sd'
        if os.path.isfile(sd): os.remove(sd)
        arcpy.StageService_server(sddraft, sd)

        print 'Publish service...'
        arcpy.UploadServiceDefinition_server(sd, connection6546543)
        print "Service successfully published"
        print 'Clean sd'
        if os.path.isfile(sd): os.remove(sd)
    else:
        print "Service publication failed: Errors during analysis."

    print arcpy.GetMessages()

    print 'Clean sddraft'
    if os.path.isfile(sddraft): os.remove(sddraft)

Detailed solution

The detailed solution I developed with the answers is the code snippet below. It has to be inserted between the draft service creation and draft service analysis. It modifies some parameters (UsageTimeout, KML off and WMS on, etc.) but can easily be adapted to others.

import xml.dom.minidom as DOM

(...)

# change service parameters
doc = DOM.parse(sddraft)
keys = doc.getElementsByTagName('Key')
for key in keys:
    if key.firstChild.data == 'UsageTimeout': key.nextSibling.firstChild.data = 6000
    if key.firstChild.data == 'WaitTimeout': key.nextSibling.firstChild.data = 60
    if key.firstChild.data == 'IdleTimeout': key.nextSibling.firstChild.data = 20000
    if key.firstChild.data == 'MinInstances': key.nextSibling.firstChild.data = 1
    if key.firstChild.data == 'MaxInstances': key.nextSibling.firstChild.data = 1
services___ = doc.getElementsByTagName('TypeName')
for service__ in services___:
    if service__.firstChild.data == 'KmlServer':
        service__.parentNode.getElementsByTagName('Enabled')[0].firstChild.data = 'false'
    if service__.firstChild.data == 'WMSServer':
        service__.parentNode.getElementsByTagName('Enabled')[0].firstChild.data = 'true'

# save changes
if os.path.exists(sddraft): os.remove(sddraft)
f = open(sddraft,"w")
doc.writexml(f)
f.close()


texts: QGIS DBManager Query between Shapefile and PostgreSQL Table

texts:
I have an instance of DB Manager running and I have TWO connections present in the left pane of the DB manager form.

First connection is to a static shapefile on my local file system
Second connection is to a PostgreSQL DB-table

I would like to run a compound query between the two in DB Manager.
However, The problem is when you instantiate an SQL window it is either associated with one or another connection. Any suggestions.
Is there a way to perform such a query with in QGIS Manager.
I am running QGIS v2.18.11 and we do not have any shp plugins, availabl e and nor the means to download them


texts: Adding Google map layer in Openlayers - what about this function?

texts:
How to add Google map as a layer (something like  "var gmap = new OpenLayers.Layer.Google("Google Streets", {visibility: false})"?) to this mymap.js function i already have:

var map = new OpenLayers.Map({
    div: 'map',
//  allOverlays: true,
    theme: null,
    projection: sphericalMercatorProj,
    displayProjection: geographicProj,
    units: 'm',
    numZoomLevels: 18,
    maxResolution: 156543.0339,
    maxExtent: new OpenLayers.Bounds(
        -20037508.34, -20037508.34, 20037508.34, 20037508.34
    ),
    controls: [
        new OpenLayers.Control.Attribution(),
        new OpenLayers.Control.Navigation(),
        new OpenLayers.Control.PanZoom(),
        new OpenLayers.Control.LayerSwitcher()

    ],
    layers: [
        new OpenLayers.Layer.OSM('OpenStreetMap', null),
        sundialsLayer,
        sprintersLayer,
        poisLayer
    ],
    center: new OpenLayers.LonLat(0, 0),
    zoom: 2
});

Thanks.

Note: suggested answer is fine -

    layers: [
        new OpenLayers.Layer.OSM('OpenStreetMap', null),
        new OpenLayers.Layer.Google(
        "Google Streets", // the default
        {numZoomLevels: 20},
        )

        sundialsLayer,
        sprintersLayer,
        poisLayer
    ],


texts: Problem with using Fiona to write out GeoJSON objects to GeoPackage layers

texts:
I'm attempting to use fiona to write out a number of GeoJSON objects to different layers in a single GeoPackage.
Here's an example json.dumps() representation of the GeoJSON object containing one feature:
{"type": "FeatureCollection",
 "features": [
    {"type": "Feature",
     "geometry": {
         "type": "Polygon",
         "coordinates": [[bunch of coord pairs]]},
     "properties": {
         "dataType": "Shapefile",
         "fname": "cb_2018_us_state_500k.shp",
         "path": "file:///C:\\Data\\Foundation\\US\\cb_2018_us_state_500k",
         "native_crs": 4269,
         "lastmod": "2019-04-15T15:55:28"}
}]}

This feature collection, call it feature_collection_a, is stashed as one of many, sharing the same schema, in a dict I'm iterating over to write the outputs:
layers = {'level_00': feature_collection_a,
          'level_01': feature_collection_b}

So I try to do (assuming I've imported from_epsg from fiona.crs):
path = 'C:/path/to/someoutput.gpkg'
driver = 'GPKG'
schema = {'geometry': 'Polygon',
          'properties': OrderedDict([
              ('dataType', 'str'),
              ('fname', 'str'),
              ('path', 'str'),
              ('native_crs', 'int'),
              ('lastmod', 'str')])}

for k, v in layers.items():
    if v:
        with fiona.open(path, 'w',
                        schema=schema,
                        driver=driver,
                        crs=from_epsg(4326),
                        layer=k) as outlayer:
            outlayer.write(v)

This gives me:
Record does not match collection schema: dict_keys([]) != ['path', 'fname', 'dataType', 'native_crs', 'lastmod']

Since I am using an OrderedDict() to construct both the properties element (in the raw geojson object represented by v) and in schema, I am confused as to why this is failing.
When I try it with geopandas it works fine:
import geopandas as gpd

for k, v in layers.items():
    if v:
        gdf = gpd.GeoDataFrame.from_features(v)
        gdf.crs = 'EPSG:4326'
        gdf.to_file(path, driver=driver, layer=k)

And since geopandas uses fiona for its r/w ops, it tells me the problem is not with fiona but with how I'm (mis)using it.
I could just use geopandas to do this, but (A) I feel like that's lazy and (B) relying on a very high level package like GeoPandas makes my script more fragile than I'd like.


texts:
In QGIS 3.4 I can run algorithms under the General setting "Invalid features filtering" -> "Ignore features with invalid geometries"
However, when I try Processing through Python functions like the qgis:fieldcalculator it just breaks.
Is there a way to make my QGIS python interpreter recognize that I do not care for invalid geometries?

Here is an example of a working code with valid geometries which breaks when invalid :

import os
import sys
import gdal
gdal.UseExceptions()

from qgis.core import QgsApplication
from qgis.analysis import QgsNativeAlgorithms

QgsApplication.setPrefixPath(os.path.join("C:", os.sep, "OSGeo4W64", "apps", "qgis"), True)
qgs = QgsApplication([], False)
qgs.initQgis()

sys.path.append(os.path.join("C:", os.sep, "OSGeo4W64", "apps", "qgis", "python", "plugins"))
import processing
from processing.core.Processing import Processing

Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

def add_area(input_vector):
    out = "pathtooutput.shp"
    params = {
                'INPUT': input_vector,
                'FIELD_NAME': 'area',
                'FIELD_TYPE': 0,
                'FIELD_LENGTH': 10,
                'FIELD_PRECISION': 3,
                'NEW_FIELD': True,
                'FORMULA': ' $area ',
                'OUTPUT': out
                }
    res = processing.run("qgis:fieldcalculator", params)
    print(res)
    return out

inp = "path_to_testvector_with_invalid_or_valid_geometries_for_test.shp"

add_area(inp)


texts: Using spatial indexes on point to polygon intersection query in SQL Server

texts:
I have a table (tab_point) with about 5 million points and another table with a few thousand polygons (tab_polyg). I need to find what points intersect with what polygon (to establish a count of points per polygon). I am really not sure how to use spatial indexes to speed up the query.
I am using this query (it's been running for hours...):
SELECT b.POLYG_NAME as POLYG_NAME, count(*) as POINT_COUNT
INTO [result_table]
FROM [tab_point] a WITH (INDEX(tab_point_spa_idx))
JOIN [tab_polyg] b
ON a.Shape.STIntersects(b.Shape) = 1
GROUP BY b.POLYG_NAME;

As we can see in the query I am using a spatial index for the point table. But I am doing this without much knowledge really.
Generally speaking, if I have a spatial query establishing a relationship between 2 tables, how can I know which one should be using the spatial index?
Should it actually be both?
And also, would it make any difference if I change the order of the STIntersect command? ("ON b.Shape.STIntersects(a.Shape) = 1" instead)


texts:
I am using QGIS 3.8 and installed python 3.9 where my Scipy is linked (I guess). The thing is my QGIS has a warning error that python that says:
Semi-Automatic Classification Plugin: Error. Check Python Scipy installation for
the Semi-Automatic Classification Plugin

When I checked the python version in QGIS it says there that its 3.6. Is there a way to change the python version in QGIS so it recognizes the Pythn 3.9 where I guess I have my Scipy fixed. I'm using semi automatic plugin to process some images.
I am not a programmer so I don't know any of this I'm so sorry.


texts: Using the same backlink and cost-distance rasters for different source/destination points over the same geographic region

texts:
I am conducting a cost distance analysis for many points on the same map (roughly 350 sets of source/destination points, creating 350 different cost distance paths in total). The map is large, at approximately 3.0*10^5 km^2. I have created the cost surface for the geographic region, as well. All that needs to be done is to run the information in the Cost Distance tool with a loop via Arcpy. This task will be time consuming, though, and I am wondering if there are any "shortcuts" I can take to reduce computation time.

Since I am using the same map to conduct my cost-distance analysis, can I somehow use the same backlink and same cost distance rasters created with the first source/destination points for the remaining ~349 source/destination points, produce incorrect results?


texts:
I'm trying to make a small change in this QGIS plugin (https://github.com/camptocamp/QGIS-SpreadSheetLayers), but I did not succeed.
I have an ods spreadsheet in which the coordinates can come in meters or in degrees, minutes and seconds format (DD MM SS,SSS), separated by spaces.
In the first case the plugin imports the spreadsheet without problems, despite the decimal separator being comma. However, in the second case the plugin adopts only the value of degrees for latitudes and longitudes. I'm debugging the plugin for two weeks in an attempt to make this change but came to fail.
Can anybody give me a hint where I can make this change?

First case

Second case


texts:
I am running a Salesman analysis and get the following error message:

Terminal at node [1860] cannot be connected to terminal at node [1923]

I have been able to identify one of the nodes by a process of elimination, but is there a way for me to identify node #1860 and #1923 easily and visually?

This relates to my other question Orphan road segments


texts: Probability of occurrence in numpy raster with GDAL

texts:
I generate a lot of different flood raster maps with a model. After each model iteration, I want to update a numpy array with the occurrence of a certain value - or  more precisely: the probability of occurence of a certain value through all model runs. With ArcPy that'd look like that:

prob = Con(value < 1, prob + float(1) / float(number_of_runs), prob)

I want to do this with GDAL however. But I am not sure how to access and iterate through every cell. Do I need to flatten the array first? My raster maps are quite large, so performance is an issue!

In the code below is the link to a sample DEM. What I want to do here: for each of the 20 model runs, update the output raster where the value is < 2350, else use the value from the output raster. As this is only a static elevation raster, the cell values don't change, so I'd expect an output raster with 1 (100%) wherever the elevation is < 2350, everywhere else it's 0 (0%). What's relevant in the code example is the part right under #The actual calculation

# get small DEM from https://github.com/GeospatialPython/Learn/blob/master/dem.zip

import numpy as np
from osgeo import gdal, osr

# Import raster & transformation to numpy array
ds = gdal.Open("dem.asc")
prj=ds.GetProjection()
scrs=osr.SpatialReference(wkt=prj)

data = np.array(ds.GetRasterBand(1).ReadAsArray())

#The actual calculation
dataOut = np.empty(shape=[ds.RasterYSize, ds.RasterXSize]) # create empty np array

#Calculate probability of occurrence
nruns = 20
for i in range(0,nruns):
    for j in data:
        if data[j] < 2350:
            dataOut[i] += 1.0 / nruns
        else:
            dataOut[i]

# Create output raster
output_raster = gdal.GetDriverByName('GTiff').Create('dem.tif', ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32)  # Open the file
geotransform = ds.GetGeoTransform()
geotransform = geotransform[0],geotransform[1],0,geotransform[3],0,geotransform[5]
output_raster.SetGeoTransform(geotransform)
band = output_raster.GetRasterBand(1)
band.SetNoDataValue(-9999)
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS(scrs.GetAttrValue('projcs'))
output_raster.SetProjection( srs.ExportToWkt() )
output_raster.GetRasterBand(1).WriteArray(dataOut)
output_raster = None


texts: Building Python add-ins without needing to restart ArcMap?

texts:
Using the Python add-in wizard I need to restart ArcMap every time I want to test something.  For example, when I've got syntax errors or something little that I need to troubleshoot.

What approach would you recommend to negate the need to start ArcMap repeatedly?

It has been noted in comments that this question has already been asked and vaguely answered; however it's totally unclear to me how you integrate this into your script.

Can someone give me an example?

I've tried to include the following script as noted in the other question linked here.

import os

esri = "C:\Program Files\Common Files\ArcGIS\bin\ESRIRegAddIn.exe"
cwd = os.getcwd()

mapdoc = "C:\Data\Shaun Pilot Area\PilotAreaMXD.mxd"

try: os.system("TASKKILL /F /IM ArcMap.exe")
except: pass

os.system(os.path.join(cwd, "makeaddin.py"))

os.system('"{0}" {1} /s'.format(esri, os.path.split(cwd)[-1] + ".esriaddin")

os.system(mapdoc)

Here's the output I get

ERROR: The process "ArcMap.exe" not found.
'C:\Data\Shaun' is not recognized as an internal or external command, operable program or batch file.
'C:\Program' is not recognized as an internal or external command, operable program or batch file.
'C:\Data\Shaun' is not recognized as an internal or external command, operable program or batch file.


texts: Is there a way to represent just the "data" footprint in a tile index?

texts:
Is there an easy way to get just the "data" footprint for a raster. I have 3 DEMs and a tile layer. I'll be using the tile layer for referencing tiles of contours and ultimately I want to flag tiles along the edge (with a field) as "incomplete" if the entire tile isn't covered by data.

I tried to use gdaltindex because I remembered using it to create reference layers for imagery and I figured it would do what I want but the result includes all the NoData in the extent of the image and I can't seem to get around that. I was hoping to do a quick selection of the tiles that are "completely contained" and then switch selection to get the ones that are outside the "data" area.

I figure I can (as a last resort) just render the DEMs as all values of 1 and mask out the NoData values and then convert to polygon but I was hoping there was an easier way?

I can make use of a QGIS or ArcMap solution if someone has one.


texts: Is there a way to convert an .lasd to multipoint? I created this .lasd from .laz files using the 'Convert LAS tool' in ArcGIS Pro (.las were unavailable).


texts: projecting onto Cartographic coordinates using Gauss Conformal Mercator in python

texts:
The shapefiles and boundaries I work with are on WGS 84 referenced to SA:

GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]]

I need to project them onto a plane to make area calculations and intersect shapes. I'm stuck with a few questions:

local official projection is defined by Gauss Conformal Mercator. How do I project files with this map definition in Python?
here are a number of definitions: http://spatialreference.org/ref/epsg/2047/, as well as http://epsg.io/?q=south%20africa.

Which package/wrapper/function do I choose to do the heavy-lifiting?

I read that planar projection in SA are done in a variety of ways- Which of the Lo/Wg do I choose for the projection to cover the entire span of the country?

The mapping needs to maintain the areas. Does tha GCM do that?
Is Gauss-Krüger the same as Gauss Conformal?


texts: Setting PYTHONPATH for compiled QGIS version on Linux

texts:
I'm trying to set the PYTHONPATH in order to be able to load the qgis classes in a python console (not the QGIS Python one).

I followed some resources I've found (PYQGIS Cookbook, and others) but I get still the same error message:

>>> import qgis.core
ImportError: No module named qgis.core

I'm pretty sure the problem is that I have a compiled version of QGIS and of course the PATHs suggested are different.

The folder of my QGIS is:

/home/user/QGIS-master

I tried to add the following to the .bashrc file in order to have the PATH but nothing changed (still the error of above):

export PYTHONPATH=/home/user/QGIS-master/python/

Someone has an idea of which is the correct PATH I have to export?


texts:
I see lots of references for Tasseled Cap transformation for Landsat-8 Top of Atmosphere reflectance.
Is the Tasseled Cap approach applicable for Landsat 5, 7, and 8 surface reflectance data as well?

Landsat 5, 7, and 8 surface reflectance Tasseled cap transformation


texts: Different level of detail(layers) on different zoom level at  OpenLayers map

texts:
I want display different detail(layers) on different zoom level on OpenLayers map. I mean when zoom in should display lot of layers and when zoom out those unsuitable layers should hide automatically. Different layers should display on different zoom level ,not using layer switcher.

Can this do using OpenLayers?

If not what is the open source solution?
If yes, pls give some example or links to get some idea

Thanks for reading and answering


texts: Importing .csv to QGIS gives CRS undefined (but not always)

texts:
I wonder why some .csv files when imported (with long and lat) can not set a proper CRS system and other do (CRS undefined). For example:

If I export my data in csv with Microsoft Excel it does not set the CRS and can not be fixed.
when I export with Google Sheets it does work.

What is the problem in this case? I am trying to export data with pandas python but can not figure out this problem? I am using QGIS, 2.12.1-Lyon, on MAC.


texts: Get data from WFS and present in using d3js - correct projection

texts:
I'm trying to get data for map of Poland. I found WFS service that provides GML files. It is located under:

http://sdi.geoportal.gov.pl/wfs_prg/wfservice.aspx?request=GetCapabilities&service=WFS

I'm able to all provinces of Poland using this request:

http://sdi.geoportal.gov.pl/wfs_prg/wfservice.aspx?request=GetFeature&version=1.1.0&service=WFS&typename=gmgml:WOJEWODZTWA

My problem is that data from this server is in EPSG:4326 projection, when d3js is using U.S.-centric (d3.geo.albersUsa).
I can convert GML to shapefile (using QGIS) and then to TopoJSON or GeoJSON (using http://www.mapshaper.org/), but I need to know what projection I must use to display my map correctly.


texts:
I am new to scripting, as the following question will indicate, but I am trying to build a QGIS python plugin. The plugin has two comboboxes, the second one should show a list of fields from the layer selected in the first combobox. Here is my script:

def run(self):
    """Run method that performs all the real work"""

    layers = self.iface.legendInterface().layers()
    layer_list = []
    self.dlg.comboBox.clear()
    for layer in layers:
        layer_list.append(layer.name())
        self.dlg.comboBox.addItems(layer_list)
    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    print result

    index=self.dlg.comboBox.currentIndex()
    layer = self.dlg.comboBox.itemData(index)
    Field_list = []
    self.dlg.comboBox_2.clear()
    for field in layer.pendingFields():
        Field_list.append(field.name())
        self.dlg.comboBox_2.addItems(layer_list)
    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    print result

But when i run this i get an error:

  Traceback (most recent call last): File
  "C:/Users/carmel.han/.qgis2/python/plugins\zoomtoselect\zoomto_select.py",
  line 201, in run layer = self.dlg.comboBox.setCurrentLayer()
  AttributeError: 'QComboBox' object has no attribute 'setCurrentLayer

I suspect i am missing a module. Anybody knows if that is in fact the problem and how do i fix it?
I have 4 modules in my PyQt4 folder: QtCore, QtGui, QtNetwork & QtXml.


texts: ArcGIS 10.2.1 arcpy handling quote mark styles

texts:
I have been asked a question that I believe is due to a mix a quote marks however I cannot specifically say why this is so in the example below. I wondered if anyone has had a similar issue.

A small amount of code has been written (with the string being copied and pasted), which ends in a invalid syntax error message. e.g.

>>> for dataset in datasetList:
...     print dataset + " " + “Dataset”
...
Parsing error SyntaxError: invalid syntax (line 2)

Above you can see there appears to be a mix of quote marks in terms of font but I cannot see why that is of importance as they are both double quotes.

Why can't the software interpret different fonts of the same keyboard symbol?

My responce to the user is to ensure all quote marks are consistent either ' or ".

To test I edited the code and it works without a syntax error...

>>> for dataset in datasetList:
...     print dataset+ " " + "Dataset"
...
Culture Dataset
Assets Dataset
Pipes Dataset
Wells Dataset
AerialPhoto Dataset

I also tested using different styles of quote marks (e.g. ') to see if the issue was due to a mix of quote marks being used. I ended up with two successful excutions of the code:

>>> for dataset in datasetList:
...     print dataset, 'Dataset'
...
Culture Dataset
Assets Dataset
Pipes Dataset
Wells Dataset
AerialPhoto Dataset

and

>>> for dataset in datasetList:
...     print dataset + ' ' + 'Dataset'
...
Culture Dataset
Assets Dataset
Pipes Dataset
Wells Dataset
AerialPhoto Dataset

Interestingly this does not seem to be a problem when using subsequent versions of ArcGIS i.e. 10.2.2 and 10.3 but maybe this is a consequence of proper quotation marks being used as supposed to software handling this specific problem.

Has anyone else come across this problem and if so what is the reason for the computer not seeing " and ” as the same symbol?