texts: How to split geometry at points using Arcobjects?

texts:
I am attempting to split the geometry of a line feature for further manipulation, without modifying the feature itself.  I am using an IPointCollection, consisting of points that intersect the line, to split the line, where I can find the midpoint of these line segments.  My code:

IFeature m_feature = some polyline feature;
IFeatureClass pointFeatureClass = some point feature class;

// storage for point geometry
IGeometry multipoint = new Multipoint() as IGeometry;
multipoint.SpatialReference = m_feature.Shape.SpatialReference;
IPointCollection intersectingPoints = (IPointCollection)multipoint;

// get intersecting points
sFilter.Geometry = (IGeometry)m_feature.Shape;
sFilter.GeometryField = pointFeatureClass.ShapeFieldName;
sFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;

IFeatureCursor pCursor = pointFeatureClass.Search(sFilter, false);
IFeature thisFeature;
object Missing = Type.Missing;
while ((thisFeature = pCursor.NextFeature()) != null)
{
    IPoint pnt = thisFeature.Shape as IPoint;
    pnt.SpatialReference = thisFeature.Shape.SpatialReference;
    intersectingPoints.AddPoint(pnt, ref Missing, ref Missing);
}

MessageBox.Show(intersectingPoints.PointCount.ToString()); // returns correct number (i.e. greater than 1)

// split m_feature polyline at intersecting points into individual polyline objects
IPolycurve2 polyCurve = (IPolycurve2)m_feature.Shape;
IEnumVertex splitPoints = intersectingPoints.EnumVertices;
IEnumSplitPoint enumSplitPoint = polyCurve.SplitAtPoints(splitPoints, true, true, -1);

// check if split occurred
if (enumSplitPoint.SplitHappened)
{
    MessageBox.Show("Split Happened") // always shows as split
}

// new geoCol for polycurve
IGeometryCollection geoColl = (IGeometryCollection)polyCurve;
MessageBox.Show(geoColl.GeometryCount.ToString()); // prints 1 as if split didn't occur

Referencing Splitting zAware polyline by points using ArcObjects with C#? and GeoNet.

What am I doing that might be causing the polycurve object to not be split?  Or is there a better way to accomplish what I am doing?


texts: How to determine which GDAL version rasterio uses?

texts:
I want to resample a georeferend raster .tif file with Python using rasterio.

I found this:
example

and here is a description of the
resampling methods

I have to use the max method which is available in GDAL version 2.0.0 and higher. I have installed GDAL version 2.4.4 via brew on my Mac OS system.

When I run the following code (i.e. almost the example code):

raster = Path('dsm_EPSG_4839_cropped.tif')

import rasterio
from rasterio.enums import Resampling

upscale_factor = 1/(2**5)

with rasterio.open(raster) as dataset:

    # resample data to target shape
    data = dataset.read(
        out_shape=(
            dataset.count,
            int(dataset.height * upscale_factor),
            int(dataset.width * upscale_factor)
        ),
        resampling=Resampling.max
    )

    # scale image transform
    transform = dataset.transform * dataset.transform.scale(
        (dataset.width / data.shape[-1]),
        (dataset.height / data.shape[-2])
    )

I am getting the following error:

Segmentation fault: 11

Is there a way to determine the GDAL version the rasterio package uses? Is there a way to change that GDAL version rasterio is using?


texts: Using PostGIS Geohash function to run through column of coordinates and dump results in new column?

texts:
I'm new to databases.

How can I use the PostGIS Geohash function to run through a column of coordinates and dump results in new column?

Here is the postgis documentation on it

Here's what my data looks like:


texts: Improve performances on distance coverage routing - PostgreSQL

texts:
I have to create layers of land coverage from one point to anywhere around on a given street network for a given distance (ex: FROM house, DISTANCE 30000m)

An important factor to understand; I need the result to be exactly superposed over my network, and the coverage to represent only the parts of the road (let say, the addresses) that can be reached from the network for the given distance. Nothing more.
I need to be precise, because it will be used for emergency services and efficiency analysis over the covered network.
And to add to the fun, the tool must not only give the right answer, but do it in a decent length of time, as it will ultimately be open to other uses on a public web interface (so no, I don't only need an application that will generate the right routing over a night of work!)

The point is not to find the shortest way to get to a known point, so basic pgrouting tools can't get me a direct answer.

I developped a working algorithm, and it works well on long routes, but as soon a I enter a "real" network (as for a city, or even a village), it crashes... well, it simply run for hours.

Basically, my algorithm goes like this:

First, we link our network table to our previous results to make sure we don't take a route that as been looked uppon, unless we have more "remaining distance" to cover at this point, which mean that our new route is better.

  FOR a in EXECUTE 'SELECT aa.ogc_fid, aa.'||col_contr||' as cost, aa.the_geom, aa.source, aa.target

   FROM '||table_name||' aa LEFT JOIN temp res ON res.ogc_fid = aa.ogc_fid AND res.dist > '||distance||'

    WHERE (aa.source = '||debut||'
            OR aa.target = '||debut||')
            AND (ARRAY[aa.ogc_fid] &&  ARRAY['||array_to_string(fromSegment,',')||']) = false
            AND res.ogc_fid is null'

      LOOP

After that (I'll short the code),
IF(distance > a.cost) THEN

We insert the informations in the result table, and launch back the function from this node, which developps a tree.

ELSE (we don't have enough "remaining distance" to cover the line)

We return a st_line_substring of the lenght we can cover, and the function stops

So it's a reccursive function that tries to cover every possible routes from the startpoint to the given distance...

Any idea about how to improve it's speed?
I think the answer would be to modify the fist FOR statement... but after hours of brain boiling, I have to admit that I can't do it alone!

Thanks

Some new details

Suppose this fictive (and horrible, sorry, not on my computer) network, with the blue point as a startpoint for a given distance.

Driving_distance would give a result such as this (red dots being end vertex for the distance, and green ones intermediate vertex on the path).

The problem I face here (but still, not on my computer, can't try it tonight, maybe I'm wrong) is that I need a very precise and stuck-on-the-network polyline result.
With these dots, I guess I could recreate the path, using some kind of point_to_line function, but it would not work for complexe polylines that are not straight from startpoint to endpoint (like, if the road makes an "S"...).

Same problem with isocrone:

The coverage here is not what I need, simply because it covers spaces that cannot be reached on the real network (like the segment in the pink ellipse, which is pass the terminal red dot).

(Sorry about the pictures, I'll probably give you some real screenshots of the problem tomorrow)


texts: Plot overlapping buslines (polylines) with leaflet in R

texts:
Long story short: I´d like to plot bus lines on my GeoJSON data with Leaflet but need to get rid of the overlap.

I´ve been skimming the forum for an answer to my question but, haven´t found any yet. I was trying to use the polylineoffset.js https://github.com/bbecquet/Leaflet.PolylineOffset for Leaflet but couldn't figure out how it works. It may be due to the fact that I´m trying to plot the map in R Studio only.

Code:

leaflet()%>%
addTiles()%>%
setView(lat = 51.545483,
        lng=9.905548,
        zoom=18)%>%
addMarkers(lng = Ein$Longitude,
           lat = Ein$Latitude,
           popup = str_c(as.character(Ein$Ein),"Einstiege:",as.character(Ein$Einstiege),
"Ausstiege:", as.character(Ein$Ausstiege),sep="
")
          ,label =  as.character(Ein$Ein))%>%
addPolylines(data=lines1,
             color=~pal(lines1@data$`(pers = n())`),
             fillOpacity = 1,
             popup =str_c(as.character(lines1@data$ref), "Durchschnittliche
Stationen", as.character(lines1@data$`(Stationen =
mean(round(StationenGesamt, 2), na.rm = TRUE))`),"Durchschnittliche
WLQ:",as.character(lines1@data$WLQ),sep="
"))%>%
addLegend(position="bottomright",
          pal=pal,
          values=lines1@data$`(pers = n())`,
          labels)%>%
registerPlugin(polylineoffset_Plugin)

How can I do this?


texts: Print elevation of different coordinates from .hgt file

texts:
I am trying iterate through the data(elevation) from .hgt file. I have attached the codes which i used to get the elevation from the file.

from gmalthgtparser import HgtParser
hgt_file='N09E076.hgt'
with HgtParser(hgt_file) as parser:
    alt = parser.get_elevation((9.666673719181665,76.30757542014408))
    print(alt)

The output of the above code will be : (1200, 1108, 7)

So here the input is given by the user with exact longitude and latitude. I want it to iterate through all the coordinates under some range and print each coordinates elevation and it should be saved as an excel file.


texts: Is it possible to publish a geoprocessing services that uses another dll?

texts:
I have a functional model (Note:  not something from ArcGIS Model Builder) that is distributed as a compiled dll that I want to access via ArcGIS.

Right now I am writing python wrappers to communicate with the dll (convert the information in the features to the correct format for the dll, then convert the model output into formats for ArcGIS and my features).

Once I have finished that and have it running locally, is it possible to deploy the python code as a geoprocessing service, and the dll on the same server as the ArcGIS server, and have it all work?


texts: converting coordinate/reference system on IDRISI Selva

texts:
i had converted my shp.coordinate/reference system from plane to UTM 48S in IDRISI Selva by change it directly at its layer properties. And after, theres message said that the changing is only for temporary and not change its origin properties.

Any suggestion ?


texts: ArcGIS Pro 2.5 - Raster stretch - change minimum maximum values

texts: ArcGIS Pro 2.5 - When I set a raster to stretch symbology and use a minimum maximum for the value range, how/where do I change the minimum and maximum values?  In ArcMap there used to be a little tick box than would let you override the minimum and maximum, where has that functionality gone? (assuming it is even in ArcGIS 'Pro')


texts: ArcScene fails to export multipatch to 3D?

texts:
I'm having a problem with the exportation in ArcScene. I want to export multipatch layer consisting of two spheres to the VRML (*.wrl). But ArcScene crushes down every time. I tried different options like JPEG texture quality etc. but it didn't work. I want to export this multipatch with a DEM. Exporting DEM alone works. But multipatch is causing some trouble.

Any ideas why?


texts:
I am trying to find a way to query my two tables in my database. So the first table contains WiFi hotspots and the other table contains city zones. So the.geom and geom are geometry columns and I am trying to see if the WiFi geometry is within the zone geometry, also I am trying to select some range within the zones. Here is my query,

 SELECT  distinct a.geom
FROM wifi a, zones b where ST_Within(a.thegeom, b.geom) and
b.zone BETWEEN '101' AND '108'
OR b.zone BETWEEN '301' AND '303'
OR b.zone = '306'
OR b.zone = '402';

After this query is done, instead of returning the geometry within the zones range, all geometries are outputted. I'm a bit annoyed and frustrated with this.


texts: How to extract values from multiple rasters and polygons for given coordinates

texts:
I have roughly 30  predictor variables in raster and two polygon shapefiles format that I am trying to extract values from for a given response variable coordinates.

I have already done the following:

reprojected all data into same projection (GDA 1994)
resampled all rasters to have same resolution

Now my specific question is as follows:

What function will allow me to extract values from all the predictor variables for given coordinates and output into a .csv file?

Just to be clear, the following is the output format I need:

X   Y   Insect.Abundance Predictor.1 Predictor.2 .... Predictor.x

My predictor variables are in rasters and polygons. What I want to do is extract a value for each X Y point provided by my response variable. I have not done anything with the polygons besides check for correct projection.


texts: Creating O-D cost matrix using ArcGIS Pro with routes from network data and not just straight lines

texts:
I want to create a O-D cost matrix which will relate the origins and destinations by a transportation network layer and will not just connect by a straight line.
Can you give me any insight on this?


texts: How to export Polygons to CSV with coordinates?

texts:
I am very new to QGIS, so please bear with me.

I have a ".shp" file, which is contains a set of boundaries for suburbs of Australia. I opened it to QGIS and I would like to export it into a CSV file so that the final file would be like following.

POSTCODE, LATITUDE, LONGITUDE

When I select "Layer > Open attribute table" I can see the "postcode" column.

I had a look at a similar question here

Getting list of coordinates for points in layer using QGIS?

Which told me to simply go to "Layer > Save as" and then type "GEOMETRY=AS_XY" or "GEOMETRY=AS_XYZ" in the "Layer" section of the "OGR creation options".

I followed exactly, but the latitude and longitude values are not being returned.


texts:
My QGIS installation has python 2.7 by default. I am not able to update this python 2.7 to an upper version of python ( e.g. python 2.7.11).
I searched the web, however, I do not want to make changes without being sure. Is it ok if I install the new python and set the PYTHONPATH environment variable as below:

set PYTHONPATH=c:\qgispath\python

Edit: I mean to update python 2.7 to 2.7.11 and not 2.11 ( Sorry for the typo). I updated the question.


texts:
I have a few (~50) PBF files came from SRTM conversion using phyghtmap.

How to merge all of them before importing with osm2pgsql?


texts:
I want to create a fishnet based on the extent of individual polygons from a list of polygons. I plan on processing each resulting fishnets further.

I am maintaining a worldwide data set on Quarter Degree Grid Cells split into country level grids. I have previously used a python script to generate the grid cell sets. My workflow has been as described in an article on my blog. The point which breaks my workflow is first generating a world level fishnet, certainly  when I am operating on a 1/8th degree cell level and below. I want to create the necessary fishnets on a country level instead. To support that workflow I need to know the following:

Does FME provide functionality to create fishnet polygons with a defined lon/lat starting point? And if so does that method accept scripted inputs so that this can be done based on a vector data source with n countries/regions defined?


texts: Convert transparent black to white with only GDAL

texts:
I have an RGBA foo.tif with NoData Value=0 on all 4 bands.
Unfortunately values that are transparent black (0 0 0 0), when converting to jpg, become black instead of no data becoming white like I would expect.
I am trying to work around it with a basic color = if R+G+B==0 then white else color for each band:
gdal_calc.py -R foo.tif --R_band=1 -G foo.tif --G_band=2 -B foo.tif --B_band=3 --outfile=result.tif \
  --calc="where(R+G+B,R,255)" --calc="G*where(R+G+B,G,255)" --calc="where(R+G+B,B,255)" --overwrite

Or this "arithmetic-only" variant
gdal_calc.py -R foo.tif --R_band=1 -G foo.tif --G_band=2 -B foo.tif --B_band=3 --outfile=result.tif \
  --calc='R*((R+G+B)!=0) + 255*((R+G+B)==0)' --calc='G*((R+G+B)!=0) + 255*((R+G+B)==0)' --calc='B*((R+G+B)!=0) + 255*((R+G+B)==0)' --overwrite

But the result is not as I expect:
Input / Output

Al lot of colors (orange...) have turned white, which is the output's nodata value. And if I add --NoDataValue=0 I get this:

Which makes me think the Blue Band 3 is somehow the only one impacted by the input "nodata" value of 0. How can I fix this? I have GDAL 3.2 without the new -hideNoData flag :-(
Or better, can I make a jpg with a white "bottom" with pure GDAL (like this answer, but without rasterio because I actually want mbtiles jpg in the end)
gdalinfo foo.tif
...
Band 1 Block=390x5 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=390x5 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=390x5 Type=Byte, ColorInterp=Blue
  NoData Value=0
Band 4 Block=390x5 Type=Byte, ColorInterp=Alpha
  NoData Value=0


texts:
This is related to a previous question How to determine projection parameters when customizing a projection
 I posted.

I am trying to quantitatively evaluate scale distortion associated with choosing different projection center, center azimuth, and scale factor values for a Hotine Oblique Mercator (HOM).

1) Is the following method a reasonable approach?

Using the same concept as the spreadsheet whuber created for evaluating Albers scale distortion, create a spreadsheet filled with Snyder’s equations for the HOM (ellipsoid formula, “alternative B”, page 74 in “Map Projections – A Working Manual”).  The user inputs the chosen ellipsoid parameters (a and e), and values for the "customized" projection parameters (lat/long of projection center, centerline azimuth, scale factor, and false easting/northing). The rest of the projection constants are then automatically calculated. The spreadsheet also contains cells for each lat/long pair (in half-degree increments, or whatever increments are desired) across the projection area. The scale factor and rectified coordinates at each lat/long point are automatically calculated when changing any of the projection parameters. Now, the scale factor can be numerically evaluated 1) by computing an overall average and range of scale distortion across the projection region, and 2) the point coordinates and their associated scale factors can easily be imported into ArcMap to create a visual picture of how the scale distortion is distributed. Obviously the results are just a sample and will vary depending on how many lat/long locations are evaluated, but does this sound like a valid methodology?  Spreadsheet looks like this:

2) I've also been using a distortion analysis tool created by Michael Braymen that calculates scale (and area and angular) distortion for any given projection using an "asterisk analysis" that approximates a Tissot's Indicatrix.

The tool's Python script can be viewed or there is a also Powerpoint available that describes the tool.  I have modified the script to create 50 meter asterisk lines (i.e., a 50 meter ellipse "radius"), instead of the default 5000m.

When I compare the results from this tool against what is produced by the spreadsheet method in #1 above, the numbers do not agree very well. For example:

Sampling approximately the same number of locations (~400) across the same projection extents yields:

  Avg Scale Error using method under #1 above = 0.9997200465 (0.027995%)
  Max "compression" scale error, Method 1 = 0.9994254755 (0.057452%)
  Max "expansion" scale error, Method 1 = 1.0006580056 (0.065801%)
  Range scale error, Method 1 = 0.0012325301 (0.123253%)

  Avg Scale Error using method under #2 above = 1.0001550206 (0.015502%)
  Max "compression" scale error, Method 2 = 0.9998956844 (0.010432%)
  Max "expansion" scale error, Method 2 = 1.0010584928 (0.105849%)
  Range scale error, Method 2 = 0.0011628084 (0.116281%)

Can anyone think of a reason the results would be so different? Can I interpret the scale factor at a point (method 1) as the scale distortion of an "infinitely small circle" at that point?

Also, I am aware of the many discussions on creating Tissot's Indicatrices, so I don't need to be pointed to those ... unless there is some vetted tool out there that accurately, quantitatively evaluates distortion for user-defined regional (i.e., not global) areas, accepts the HOM, is easily implemented, and is nearly free :)?  Actually, assuming the tool used in Method 2 is accurate, it work's great for my purposes. The drawback is it takes about 9 hours to run on my PC for each evaluation.


texts:
Can someone break this script in to a few lines so that I know what I am doing?

I understand the code about finding the length of a Polyline, but how to break this into step by step?

Here is the long line script:

a[f[0] for f in arcpy.da.SearchCursor("polyline.shp", 'SHAPE@LENGTH')][0]

and as a beginner I want it something like this:

for f in cursor:
     print("The length of Polyline is", "SHAPE@LENGTH")[0]

But I am not sure if I can get a result from this, and what does [0] do here?


texts: How to do versioning while editing a WFST service published using geoserver

texts: I'm trying to edit a WFST service and need to keep different versions for quality checking and acceptance. when we edit a WFST service , it will directly change the db.could anyone please suggest a method ?


texts: Trigger function PostgreSQL/PostGIS points within polygon

texts:
I wan't to update a column in PostgreSQL that counts the number of points inside a polygon using a trigger function. The function should be fired on each insert or update of the polygons or the points layer.

CREATE OR REPLACE FUNCTION public.upd_points_in_polygons()
  RETURNS trigger AS
$BODY$
BEGIN
SELECT count(p.gid) INTO NEW.num_points FROM polygons po, points p WHERE st_intersects(NEW.the_geom, p.the_geom);
RETURN NEW;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;

Here the trigger from the table polygons.

CREATE TRIGGER t_upd_points_in_polygons
  AFTER INSERT OR UPDATE
  ON public.polygons
  FOR EACH ROW
  EXECUTE PROCEDURE public.upd_points_in_polygons();

The column num_points isn't updated, but I don't see any errors?


texts: Plotting species distribution on map in R of Europe and Africa

texts:
install.packages("rgbif")
library("rgbif")

Species.Name <- "The name of the species"

#Search GBIF and download observations
Species.key <- name_backbone(name=Species.Name)$speciesKey
locations <- occ_search(taxonKey=Species.key,hasCoordinate = TRUE, return='data', limit=20000)

#Remove localities with uncertain coordinates - meaning the observer entered approximate latitude and longitude data. This is common with very old specimens/records that were collected prior to handheld GPS devices being invented
locations <- subset(locations, coordinateUncertaintyInMeters <= 5000)
head(locations)

### **Create a Map**

#### Install the maps package to download the political boundaries for the region of interest so we can map out the data.
```{r}
library(ggplot2)
#install.packages("maps")
# Download Country Boundaries
world <- map_data("world")

Not sure how to get the entire scope of Europe and Africa show these localities on the map.


texts: How to serve raster mosaic to AutoCad Map 3D from ArcGIS Server without its Image extension?

texts:
We are looking for a way to serve all the raster datasets we've collected over the years as a mosaic we have built in ArcGIS so that the AutoCAD department can simply import the mosaic and it will show the raster data in the location they are working. We have ArcGIS Server and ArcSDE.

Has anyone accomplish this?

We would prefer to do this without purchasing the ArcGIS Image Service extension.


texts: Dealing with bridges over roads in PostGIS topology?

texts:
first time posting please ask for more details if not enough provided.

I am using OSM road data and have created postgres topology for them. I am now cleaning up the data before exporting it for network analysis.

I have encountered a problem that I am struggling to resolve for (the many) cases where a bridge goes over a road. In my topology a node has been created where they cross, but it is of course impossible to turn off a bridge on to the road underneath, so my network analysis will create incorrect routes.

In this example the node id is labelled in red and the edge id is labelled in white. For example on node 29, I need to prevent the possibility of routing along edge 43 -> edge 37. Also see node 22, this node must be preserved as edges 35 and 29 merge into edge 25, but it needs to be removed from edges 37 and 26 and replaced with one disconnected edge.

The roads have a 'layer' attribute associated to them (from OSM I believe) which ranges from -5 to 5, which might be used as a flag however I've struggled with it so far.

Also, each edge has a weight attribute that must be preserved in the final output. All I require as output is a CSV of n1, n2, weight.

The dataset is large (one country), and I'd prefer a solution in postgres, however I could explore other solutions in python, R etc.

It strikes me that this must be a common requirement when cleaning road data for network analysis but I've struggled to find any documented examples.


texts: Cannot select by attribute for new attributes added by join in ArcGIS Pro?

texts:
I'm new to ArcGIS Pro and am following some tutorials. I did a join to add some data to the attribute table of one of my layers. Now, when I go to select by attribute and look for those attributes they aren't there. You can see in the picture that the fields to the right of COMID are missing from the select by attributes pop-down window.

How do I get them to appear there?


texts: Origin Destination Cost Matrix not giving results?

texts:
I am using ArcGIS Pro. I am using a network dataset that I built. Every time I run the analysis, total_time and total_length fields render no results.
Here's an image:

What am I doing wrong?


texts:
I have this DEM file downloaded from the internet. It's a text-formatted file and it looks like this:

1440   720    30 -9999   450
-9999  -9999  -9999  -9999  -9999  -9999  -9999  -9999  -9999  -9999  -9999
-9999  -9999  -9999  -9999  -9999  -9999  -9999  -9999  -9999  -9999  -9999
...

and so on...

How do I open it in ArcMap to see it visually? Or should I open it with another software?


texts: Working with complex JSON attributes in QGIS

texts:
I am a developer and have a tool that generates a geojson asset for our data team to work with using the QGIS app. Due to the original source of the data, these features have some attributes that are represented as complex hierarchical data, impossible to project into a flat schema.

We need to be able to view / edit the attribute field as effectively a JSON blob, instead of a string literal.

What facilities are available to do this in QGIS? Are there any field editing plugins that would support such a job? Ideally some kind of popup JSON editor with document format validation and collapsible nodes, similar to how Visual Studio (Code) editors work.

Alternatively some other workflow that would enable us to manage this data in an external tool and re-inject it into the source feature.

UPDATE: Got Postgres spun up and created a JSON column from an existing text column containing JSON. QGIS hides it in the attribute table. Any further recommendations on working with JSON attributes directly in the DB?


texts: OLS tool in ArcMap 10.4: Spatially Autocorrelated Explanatory Variable

texts:
I am performing ordinary least squares regression analysis on two variables: K concentration in soil as explanatory variable, and thickness of stalk of plants as the dependent. The K-concentration dataset shows significant spatial autocorrelation (p <0.01).  The OLS results in ArcMap 10.4  return a very low R2, but the residuals are not spatially autocorrelated. The F-Statistic indicates that the model is not significant.
Could the spatial autocorrelation (clustering) in the input variable be affecting the functionality of this model?

ms


texts: Running Python script on opening QGIS project?

texts:
I have a .csv file containing financial data which is generated on a regular basis outside QGIS. Then I want to create a PyQGIS script that will add the .csv file as a layer in my project. (The idea here is I want to have the up-to-date financial data as attributes data for one of the shapefiles in the project).

Is there a way to run the script automatically as soon as I open the QGIS project?


texts: Increasing size of legend symbols in QGIS3.12.2 print layout

texts:
I have pixels with a fixed size on the map, as they have a certain resolution. Now in the legend of the print layout the symbols have the same size as the pixels on the map.
I want to increase the size of the pixels in the layout to create a legend, see the image below how it looks like now:

And what I want, and what I have with other images is a legend that looks like this:

How do I go about doing this? I tried to increase the size of the pixels using the item properties, but none of these actually increase the size of the symbol.


texts: Colorizing raster in QGIS 2.2 without using OSGeo4W command line?

texts:
Total beginner in QGIS. Following this guide right now:

https://www.e-education.psu.edu/geog585/node/695

At step 4 I created the simple txt file which I now need to use to colorize the raster (slope.tif) I created before.
Now it tells me to "open your OSGeo4W command line window" to execute some commands to do the colorizing. Where can I find this command window or a command line window for the GDAL Plugins? Do I have to install something separately? I wouldn´t mind to using command line windows.
Or is there a way to do what is done in the tutorial just by using standalone QGIS Installation 2.2 (on Win 7 64bit system)? Perhaps by using the properties menu of the raster file? If yes please give some tipps for steps 4 to 8 of the guide. I guess the rest will be no problem.


texts: Appropriate method for transforming log-transformed kriging predictions back to original scale (R)

texts:
I am attempting to run some interpolations using kriging on transformed data using the transformation log (x + 1) - because my non-negative data has many zero values. The problem I have is with back-transforming the kriging estimates back to the original scale.

I am already aware of this related question Backtransformation of kriging predictions and variances but the accepted answer does not produce satisfactory results for my data.

Here's what I am doing (the data can be downloaded from : https://www.dropbox.com/sh/xnwp3zz5abnilyo/AABRVJZ0kTmWk0T9Fcp4-bVSa?dl=0/)

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

hs1<- readOGR (".", "Hollicombe_S1_L1-5_A1.2")

#There are some negative values for the PerArIn metric which must be wrong
#for now let's assume these are 0

# replaces any negative values with 0 and remove duplicate loactions

hs1$PerArIn[hs1$PerArIn < 0] <- 0;

hs1 <- hs1[which(!duplicated(hs1@coords)), ]

#Perform autoKrige with transformed PerArIn data
PerIn.krg <- autoKrige(log1p(PerArIn)~1, hs1)

So now I want to back-transform the krige_output for PerIn.krg and this is where my issue is.

In the answer to the question referenced above the recommended procedure is to use the correction of Laurent (1963), as follows:

#estimate back to original scale following Laurent (1963)

PerIn.krg$krige_output$var1.pred <- exp(PerIn.krg$krige_output$var1.pred +0.5 *
                                          (PerIn.krg$krige_output$var1.var))

PerIn.krg$krige_output$var1.var <- exp(2*PerIn.krg$krige_output$var1.pred +       PerIn.krg$krige_output$var1.var)*
  (exp(PerIn.krg$krige_output$var1.var)-1)

But this back-transformation shows predictions well above 100 - the data is percentage data so the original scale does not extend beyond 100!

For example:

    > print(PerIn.krg$krige_output$var1.pred)
...
[297] 236.275600 231.036118 223.001912 212.319480 199.312973 184.456154
...

So I tried just using exp() as follows:

PerIn.krg$krige_output$var1.pred <- exp(PerIn.krg$krige_output$var1.pred)

PerIn.krg$krige_output$var1.var <- exp(PerIn.krg$krige_output$var1.var)

Which results in more realistic predictions, but still with some values over 100:

> print(PerIn.krg$krige_output$var1.pred)
 [497] 97.465657 100.742912 103.146009 104.504026 104.671525 103.543346

Are there any insights or suggested alternative methods (either for different transformations or back-transformation procedures)?


texts: How to find the minimum band for several Features i

texts:
I have a dataset extracted from the MODIS daily snow cover dataset, combined with a DEM. It consists of seven bands with different height classes (name of the bands). These bands consist of ones and zeros, and I needed find the band with the minimum amount of ones, and save the height (name) of that band. I need to do this for several grid cells and several days. Eventually I would like to have a table with days on one axis, and grid cel on the other, and within it the height (minimum band) at that day. I am currently trying to find the minimum band but have not managed to do so.
The code can be found here.

Is there a way to extract the name of the band with the minimum value, and save that in some data structure, for a certain day and grid cell?

Snippet of the code:

var snow = function(img){
  var e500 = img.expression('E + S',{'E':img.select(1),'S':img.select(0)}).remap([0,1,2],[0,1,0]).rename('0500');
  var e1000 = img.expression('E + S',{'E':img.select(2),'S':img.select(0)}).remap([0,1,2],[0,1,0]).rename('1000');
  var e1500 = img.expression('E + S',{'E':img.select(3),'S':img.select(0)}).remap([0,1,2],[0,1,0]).rename('1500');
  var e2000 = img.expression('E + S',{'E':img.select(4),'S':img.select(0)}).remap([0,1,2],[0,1,0]).rename('2000');
  var e2500 = img.expression('E + S',{'E':img.select(5),'S':img.select(0)}).remap([0,1,2],[0,1,0]).rename('2500');
  var e3000 = img.expression('E + S',{'E':img.select(6),'S':img.select(0)}).remap([0,1,2],[0,1,0]).rename('3000');
  var e3500 = img.expression('E + S',{'E':img.select(7),'S':img.select(0)}).remap([0,1,2],[0,1,0]).rename('3500');
return ee.Image.cat([e500,e1000,e1500,e2000,e2500,e3000,e3500])};
var B = A.map(snow);

// next a function is needed to sum and find the minimum

var unweighted = function(img){var uni = img.reduceRegions({
    reducer: ee.Reducer.sum().unweighted(),
    collection: grid,scale:500});
    var sum = ee.Feature(uni);
    return sum;
};

var sumA =B.map(unweighted);
var sumB = sumA.reduceColumns(ee.Reducer.min(7),['0500','1000','1500','2000','2500','3000','3500']);
print(sumB);


texts: Memory and Processing Efficient Approach to Read a Timeseries in Python

texts:
I am trying to store a timeseries of high resolution raster images containing NDVI values into a dictionary for calculations. I am curious as to which approach may the most efficient when loading multiple large files into memory. I currently have tried two options that are both rather inefficient.

I would loop through the folder containing all of the images and use GDAL to read the values of the images into a 3D Numpy array.

img_stack = np.empty((image_width, image_length, len(image_files)), np.dtype('f'))

for i, fname in enumerate(image_files):
    img = gdal.Open(os.path.join(fn + ("\{0}".format(fname)))).ReadAsArray()
    img_stack[:, :, i] = img
    img = None

After being loaded into the array, a loop would go to every single pixel location and retrieve the value of that location in each image.
I would do something similar to the previous approach, but instead of loading all of the images into img_stack first, I would go through the images one at a time to reduce memory usage, but it would result in a longer processing time.

Would there be a more efficient way to do this procedure, using Numpy or GDAL?

The structure of the dictionary is as follows:

z_stack = {x1:
          {y1: [z11a, z11b, z11c],
           y2: [z12a, z12b, z12c],
           y3: [z13a, z13b, z13c]},
       x2:
          {y1: [z21a, z21b, z21c],
           y2: [z22a, z22b, z22c],
           y3: [z23a, z23b, z23c]},
       x3:
          {y1: [z31a, z31b, z31c],
           y2: [z32a, z32b, z32c],
           y3: [z33a, z33b, z33c]}}

and the iteration to populate the dictionary is:

z_stack = {}
for x in width:
    z_stack[x] = {}
    for y in length:
        z_stack[x][y] = []
        for z in height:
            z_stack[x][y].append(z)


texts:
I have a csv that I am importing with lat long values amongst others. I select the correct fields for both values and ensure that both the layer and project CRS are set to WSG84(UTM30N)EPSG32630. My points should appear in the UK but are appearing at 1824,645 (south of Liberia).

I have checked my points direct in google (51.069946333,-1.701573167) and they are spot on, when I hover my cursor over the location it reads wrong. I guess I have the wrong reference system set but I understood this was the one needed for GPS.

I am using quickmapservice and googlemaps, QGIS 2.16


texts: Creating rotated local grid/fishnet in QGIS with known AMG origin & rotation angle?

texts:
I have a local grid with known origin that is in AMG (metres) and a rotation angle.

How do I establish that grid so that it can be projected into known CRS eg GDA94/ MGA Zone projection?

I located an plugin called ARKgrid that seemed like it would do the job, but it does not run in QGIS3 maderia.

Any suggestions ?

Am new to QGIS and at this stage a basic user.


texts: Reading all shapefiles in a folder using geopandas and then clip them all iteratively

texts:
I have a folder full of shapefiles, I would like to read them into geopandas and then clip them all to one bounding box iteratively, however I'm stuck on the first part which is reading each shapefile into geopandas as its own geodataframe, could anyone advise? If there is another library besides geopandas that can complete this task then that would be great also.
import glob
import geopandas as gpd
for filepath in glob.iglob('E:/folder/*.shp'):
    gpd.read_file(filepath)


texts: Calculating list of existing polygons with distances from newly created QGIS/PostGIS

texts:
I have two polygon tables (SRID=4326) in postgres: t_areas and t_sites. Once I insert through QGIS 3,8 a new area in t_areas I'd like PostGIS (by a trigger) to list all sites within a certain km range from new area and insert this list to a memo field of the t_area table, for ex.: sitename01 -> 0,3 km, sitename02 -> 0,12 km... etc. The distance to be calculated between borders, not centroids.

I've managed to compose a query on existing data like:

 select array_to_string (array ( SELECT t_sites.shortname || ' -> ' ||
 round(st_distance(geography (t_areas.geom), geography
 (t_sites.geom)))/1000 || ' km' AS distance FROM t_areas, t_sites WHERE
 t_areas.id_area=133 and round(st_distance(geography (t_sites.geom),
 geography (t_areas.geom)))/1000 < 5 ORDER BY distance ASC), ', ')

and even got result like this:

 Ванино Поле -> 0.661 km, Голубково 1 -> 3.636 km, Голубково 2 -> 4.03
 km, Голубково 3 -> 3.661 km, Задубье 2 -> 0.005 km, Заорешье 1 ->
 4.582 km, Заорешье 2 -> 4.631 km, Калгановка -> 4.385 km, Кут 2 -> 1.526 km, Наволок 1 -> 3.359 km, Петровская Горка 1 - and so on...

but I can't realize how to configure a function and a trigger to force PostGIS then automatically update the t_areas.memo field since a new area serial id_area hasn't been generated yet and so I can't point the function to that id (unlike in the query above).


texts:
Is it possible to have an orthographic projection with OpenLayers? I didn't find anything specific on the matter.

Thanks Mike, the solution is what I needed. Two questions: as I am using OL 4.6.5, I had to change the last part of the code, otherwise it doesn't work. I modified the creation of the circle in

var circle = new ol.geom.Circle([0, 0], 6370997);
var style = new ol.style.Fill({
  color: "black"
});

I imagine that the function to trigger the circle, should be changed in

map.on("postrender", function(e){
....
});

However I am not sure how to change the function internal code.


texts: How to create dynamic (depend on scale) coordinates grid in ArcGIS for Desktop?

texts:
I am working in ArcGIS Desktop. I want to dynamicly change my coordinate grid in layout view depending on map's scale. How to do that?

Currently I have to change this manually every time I zoom in or zoom out.


texts: Split a line from point but with a certain distance

texts:
I am using ArcGIS 10.2. I need some suggestion about this problem, maybe there is a certain tool, but I think no. I think some python code... but I do not know from where to start  :(

I have two shapefile layers: some polylines with points intersected. So, what I want to do automatically is to SPLIT those lines from a fixed measure (for example 10 meters) before and after this Point placed on the polyline.

The Points are on the line.

For example, if I have a 100 meters length polyline, and the point is situated in the middle, I want it split on the 40th and 60th meter.


texts: Exporting features in OpenLayers 2 to KML?

texts:
I'm new to OpenLayers (2.13). I saw in the example section that it is possible to draw features like lines, points and polygons.

http://dev.openlayers.org/examples/draw-feature.html

Quite interesting. But there is no export feature with it. I mean drawing features on the map is useless if I can't export it in KML or GML or GeoJSON.

Is there any example about exporting features in KML?


texts: Sorting table on multiple fields using ArcPy?

texts:
I'm trying to sort a table based on three fields with ArcPy (which should be only one line of code), and I had a "tool is not license" error. I just saw that sorting with multiple fields as parameters is an ArcInfo feature only. I have ArcView (soon ArcEditor).

Any ideas how I could code this to bypass that limitation with ArcPy?


texts: Obtaining same resolution of pre and post UAV-Imagery In COSI-Corr?

texts:
I am currently trying to perform a correlation between a pre UAV derived image of a  landslide with a post image to analyse displacement and direction.

The images were processed through SfM processing software Pix4D. I keep getting the following error:

Any advice?

I, unfortunately, am new to this software and have little time to grasp it.


texts:
FEMA has an updated flood WMS that doesn't seem to be working.

The URLs they say to use are:

Use the following URL in your GIS Application:
https://hazards.fema.gov/gis/nfhl/services/public/NFHLWMS/MapServer/WMSServer

The WMS Capabilities file is available here:
https://hazards.fema.gov/gis/nfhl/services/public/NFHLWMS/MapServer/WMSServer?request=GetCapabilities&service=WMS&

Here is the link to the FEMA website:

https://hazards.fema.gov/femaportal/wps/portal/NFHLWMS

In the Create a new WMS connection box I put the URL in, click ok and click connect but get the following warning message: "Failed to download capabilities: Download of capabilities failed: SSL handshake failed." I have a feeling this is an issue with authentication so I sent a message to the owner of the data and they said that their service is only supported by TLS 1.2 Standard and that I should

reach out to QGIS's tech support an inquire if their product supports TLS 1.2 and if so, how to configure you client to read the protocol.

I can't connect to the server at all. I put the URL in the URL box, click ok and then connect and get a warning message saying "Failed to download capabilities: Download of capabilities failed: SSL handshake failed" –

I am running QGIS 2.18.7 on Mac os


texts: Updating multiple attribute fields in metadata XML with ArcPy?

texts:
Updating multiple field attributes through Python in an XML metadata doc export. There are multiple fields, so the XML tags for each unique field are all the same:

      <attr>
        <attrlabl>STATE_NAME</attrlabl>
        <attrdef>State Name</attrdef>
        <attrdefs>HRSA GDW</attrdefs>
        <attrdomv>
          <udom>2010</udom>
        </attrdomv>
      </attr>
      <attr>
        <attrlabl>County</attrlabl>
        <attrdef>County Name</attrdef>
        <attrdefs>HRSA GDW</attrdefs>
        <attrdomv>
          <udom>2010</udom>
        </attrdomv>
      </attr>
      <attr>
        <attrlabl>STCOFIPS</attrlabl>
        <attrdef>State/County FIPS ID</attrdef>
        <attrdefs>HRSA GDW</attrdefs>
        <attrdomv>
          <udom>2010</udom>
        </attrdomv>
      </attr>

How do I call each unique field out in a python script to update the metadata in an XML. This is where I am:

import arcpy, sys
from xml.etree.ElementTree import ElementTree
from xml.etree.ElementTree import Element, SubElement
from arcpy import env
env.overwriteOutput = True
fcpath = r""
translatorpath = r"...\Metadata\Translator\ARCGIS2FGDC.xml"
xmlfile = r"...\FullMetaFC.xml"
arcpy.ExportMetadata_conversion(fcpath, translatorpath, xmlfile)

tree = ElementTree()
tree.parse(xmlfile)

ele = tree.find("eainfo/detailed")
for e in ele.iter("attrlabl"):
    # loop through children pulling out nodes which have a attrlabl tag
    print e.text
    # do some sort if IF THEN logic check and update as necessary
    if e.text == "STATE_NAME":
        ele.iter("attrdefs") = "NAME OF THE STATE"

tree.write(xmlfile)


texts: How do we remove labels bound to markers in google maps api

texts:
I am using:
>

ruler1 = new google.maps.Marker({
    position: map.getCenter() ,
    map: map,
    draggable: true
  });

to create a marker and then create and bind a label to it as below:

    ruler1label = new Label({ map: map });
    ruler1label.bindTo('position', ruler1, 'position');

However, when removing the marker by:
ruler1.setMap(null) which removes the marker but leaves the label behind.
I tried unbinding the label to remove it by ruler1label.setLabel(null); but it does not work.
How can I remove these labels?
I have appended a screenshot to make it clear.


texts: How do I "flood" a terrain in Cesium to represent a high water level

texts: I am thinking about using a plane object to represent a water level by rendering it as a certain height above ground using Cesium. Is there any example code that


texts: Redistribution of polygons: supplementing and filling the specific area with polygons with specific size

texts:
I'm studying the impact of green wedges (red outline on the figure) on weather parameters in Copenhagen area. What I want to do is check how spatial redistribution of future afforested areas enacted by government (green and yellow dashed polygons) would alter the weather parameters if instead of being scattered around, they would form solid, continuous green wedges.

What I'm looking for are tools/plug-ins that would allow me to redistribute the overall area of enacted forest patches and feed them into the wedges, by filling the 'empty' areas between forest (green) and urban tissue (pink), based on certain criteria (mainly wind direction and area).

Let's say I would like to supplement the most southern wedge with 5 000 h of new forest. Continuously filling the area from the inner city - outwards. Sort of like pouring specific amount of water into the glass.

Is that possible? Else all I can think of is doing it manually and scooping out polygon piece by piece to get the desired areal.

EDIT: I don't want to keep the original shapes of redistributed areas. I just need to redistribute their total area to create a continuous forest tissue.
Enacted afforestation patches account for 18477,25 h in total. This amount needs to be redistributed and it's constant.

EDIT 2 (dof1985 questions):
None of the old forest is going to be removed. The green wedges I referred to (the red outlined areas) are part of the urban development plan for Copenhagen. However, due to many reasons these wedges are disappearing and being built-up.

So project is more about supplementing/expanding the still existing forest tissue and showing its influence on weather conditions. And location is fixed - inside the (roughly) red-outlined wedges.
Criteria are a mixture of how much forest wedges already contain and how does it correlate with the dominating wind directions (e.g. dominating winds come from west and west-south but these wedges contain little greenery) and then adding new areas accordingly to get the best benefit. We are running this through a chemistry-meteorological model called HIRMAL. It's one of many scenarios of land-use change so we are keeping it rather simple.


texts: How to fix "The selected file already seems to have a world file!" when image hasn't been georeferenced before?

texts:
I am georeferencing historical aerials from 1988 in QGIS 2.2.0.
Sometimes it works perfectly and other times after I have places my points and start the georeferencer a window pops up saying:

"The selected file already seems to have a world file! Do you want to replace it with the new world file?" Yes/No.

I click yes thinking it will now take me to the transformation setting window to save my image, but nothing happens. This happens for images I haven't georeferenced before.

Please let me know what I am doing wrong.


texts:
I would like to snap these two polygons so that my city boundary lines up with the fire district boundary as much as possible.

When using the Snap geoprocessing tool, what is the best way to determine a value for the distance parameter?

City boundaries are green and fire districts are red.


texts: Adding an external WMS service using REST service with Geoserver

texts:
I'm trying to add an External WMS layer to geoserver with a PHP page, using the REST service of Geoserver. I've already read the tutorials for creating the geoserver features (http://docs.geoserver.org/stable/en/user/restconfig/rest-config-api.html), but i couldn't find the way to create a new store and layer with an external WMS service.

in other words: I have the URL of the external WMS service (=store), the layer i want to publish (=layer) and the epgs i want to use.

Is there a way to load all of these information in the geoserver schema?


texts: Using OpenStreetMap for simulation of commuter routes?

texts:
I'm seeking some advice on how to build a simulation based on OpenStreetMap Data for my masterthesis. My knowledge of spatial data is very limited, since I come from another field of science, namely information systems. I have experience with object oriented programming, as well as mySQL as relational database.

The research goal is to evaluate if applications which inform the user of the current gasoline price actually have the potential to save money. So the idea is to simulate the routes people travel between their home and work by car and when the car needs gasoline to find the nearest gasoline station and refill. As Database I will use data from the www.zensus2011.de which can tell me from which administrative district people leave and to which administrative district they travel.

Within OSM I have the following Tags that can identify the administrative regions/districts which I need

de:amtlicher_gemeindeschluessel => integer
de:regionalschluessel => integer
admin_level => 6
boundary => administrative

So my idea is to use this data to generate random routes with a certain length based on other theses which describe how far people travel between home and work. One of the main problems which I could identify from reading the OSM wiki and other online posts is that the OSM data per se is not routable and needs to be transformed for that purpose. I for that I found two programs osm2pgrouting and OSM2PO. With the latter one I managed to import the map of Germany into a PostgreSQL Database. I also imported the OSM data using osm2pgsql with the "default.style" settings into another PostgreSQL Database. Now I'm a little overwhelmed by all the data and the way how to use this data.

But let's continue with the idea. What I need to understand now is how the map data with all the information and the routing data can be connected to generate the routes.
What I thought of was: I use the administrative boundary to get a list of streets within this boundary, then randomly choose one as a starting point and do the same for the destination, calculate the route and save the route to a database for further use in the simulation. Still I need to find out a way to determine the start and destination points better, but for the beginning and understanding it would be fine with the above.

Some final words on the implementation of the simulation. Currently my plan is to use python and code the whole simulation, which I hope is feasible. I found two related works which I'm currently investigating namely MATSim and www.roaf.de. Also I'm still browsing other questions here which might help me understand the topic.

To summarize my questions:

Is the way I choose feasible?
How would SQL statements look like to connect the two Data sets?
Is it possible to save generated routes to a database?


texts: How to update multiple computers with QGIS on a network?

texts:
We have ca. 30 Computers with QGIS 2.14 and I want to update them to QGIS 2.16
for the first step and in the future simply update, when it's necessary.

The operating system is windows 7. I tried to make a configured installation like it is described in the following article http://planet.qgis.org/planet/tag/windows/
But the changed bat file doesn't work.

Has anyone experience with this?


texts: Copying fields from joined table to another table?

texts:
I'm trying to write a standalone script to copy the values from 6 fields in a feature class to another feature class based on a spatial join and name query. I am stuck at updating the values from the join table to my original table.

print "Importing modules..."
import sys, arcpy, traceback

arcpy.env.workspace = r"C:\LVS\PECO\PECO\Python\Python\Script\test_code.gdb\GAS"
arcpy.env.overwriteOutput = True

# Local variables:
GFR = "GFR"
Parcel = "PECO_TaxParcels"
GFR_SpatialJoin = "C:\\Users\\aadeyemo\\Documents\\ArcGIS\\Default.gdb\\GFR_SpatialJoin"

arcpy.MakeFeatureLayer_management("GFR", GFR)
arcpy.MakeFeatureLayer_management("PECO_TaxParcels", Parcel)
# Process: Select Layer By Attribute
arcpy.SelectLayerByAttribute_management(GFR, "NEW_SELECTION", "STREET_NAME

IS NULL")

# Process: Spatial Join

arcpy.SpatialJoin_analysis(GFR, Parcel, GFR_SpatialJoin, "JOIN_ONE_TO_ONE", "KEEP_ALL")

path_dict1 = {}

SC1 = arcpy.da.SearchCursor (GFR_SpatialJoin, ["TARGET_FID", "STHSNUM", "STSTNAME", "STSUFFIX", "STCITY", "STSTATE", "STZIP"])
for row in SC1:
   path_dict1[ row[0] ] = row[1]
   print SC1[1]
   print path_dict1
del row

ucursor = arcpy.da.UpdateCursor(GFR, ["OBJECTID", "STREET_NUMBER", "STREET_NAME", "STREET_SUFFIX", "CITY", "STATE", "ZIP_CODE"])
for urow in ucursor:
   if path_dict1.has_key( urow[0]):
   urow[1] = path_dict1[0]

This currently returns a value of zero when i try to copy values.


texts: Equivalent method for .withcolumn() for geodataframe

texts:
I am trying to use a user-defined function in pyspark in order to check whether the points in CANmep_GDF lie within the polygons in gda_GDF.

from geopandas import GeoDataFrame

def check_exist(geom):
    gda_GDF = GeoDataFrame.from_file("/dbfs/FileStore/tables/gda_000a11a_e.shp")
    exist = gda_GDF.contains(geom)
    return exist

check_exist_udf = udf(check_exist,BooleanType())

CANmep_GDF_exist = CANmep_GDF.withColumn('EXIST',check_exist_udf('geometry'))

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

AttributeError: 'GeoDataFrame' object has no attribute 'withColumn'

I think geopandas does not support .withColumn(). What would be an alternate way to write the same?


texts: Issue with Projections from WGS1984 to GCS_Amersfoort

texts:
I am using ArcGIS 10.6.1. I have data from TOP10NL that displays the road network shapefile, and I have the busline shapefile from openOV. The TOP10NL data is in GCS_Amersfoort, the coordinate system for the Netherlands. The openOV data is in WGS1984. I tried projecting the data from WGS1984 to GCS_Amersfoort (1st time) and then from WGS1984 to RD_New (the projected coordinate system for the Netherlands). I left the geograhpic transformation box blank both times, as normally it autoselects. There were no drop-down options for the geographic transformation.

The data is visualizing, but the bus line shapefile from openOV does not line up with the road network shapefile from TOP10NL as it should. I don't know what else could be wrong, but this seems to me to be an issue stemming from the bus line data being in WGS1984.

In image: Road network in GCS_Amersfoort shown in green, bus line network re-projected into GCS_Amersfoort shown in purple


texts: Carto https basemap 500 Domain not found

texts:
I have an issue with calling this tile server:

https://cartodb-basemaps-abcd.global.ssl.fastly.net/light_all/13/4505/2407.png

it says 500 (Domain not found), yet I can see the tiles on the screen. The issue I have is that I want to use the tiles in the leaflet-offline-map API, but this error causes issues and I can't download the tiles.

Any ideas?


texts: Finding calculate value or collect values tools?

texts:
I am using ArcGIS Desktop 10.5.1 and I need to use the calculate value tool or collect values tool to in ModelBuilder, but I cannot find those tools.

Do I need to download them from somewhere else?


texts: Getting business listing data for mapping for New York?

texts:
How should I go about efficiently compiling the data to map the locations of all the businesses within a given city or town? Are there any free online resources that offer searching businesses by geography and downloading their listings as tabular data which could be geocoded?

At a minimum I would just need the business name and address. The project would be in New York and the business data needs to be mapped in ArcMap.


texts:
I'm attempting to interpolate the NU national grid square of Ordnance Survey's terrain data (10m intervals). If I interpolate with IDW, I encounter no issues and the raster pixels are square in shape as they should be. When interpolating with TIN (much better for conducting viewsheds over) the NU national grid square returns non-square pixels after the interpolation process. When using the viewshed plugin this creates obscure viewshed as shown here. Can anybody else confirm you receive non-square pixels?

When trying to fix this, I was given an error "wkb access out of bounds".

This error has been appearing more and more ever since I first encountered it to the point where its occurring consistently with every TIN interpolation. I have tried the following:

restarting QGIS and the computer

different national grid squares

meticulously gone over CRS selections

checked the validity of the contours

clipped the contours to smaller areas

downgraded my version of QGIS - all at no avail.

I'm using QGIS 2.18.10. I'm conducting the TIN interpolation by navigating to Raster > interpolate and using the viewshed plugin (or would have been if the interpolation worked).


texts: Polygon layer symbology not displaying correctly in QGIS 3

texts:
I have a polygon layer with multiple attributes. I want to display a color ramp for one of those fields (which is a "real" number field calculated using 2 of the other fields).

The problem is, when I set the colors (which display fine in the "symbology tab", see first picture), they do not appear on my layer's rendering. It seems that the color for the highest value is displayed for every polygon instead of the chosen color ramp.

The display works fine for any other field, except this one (second picture).

Has something like that occured for anyone?


texts:
When a popup is opened with Leaflet a nice fade in animation is played. When you close the popup however it abruptly disappears. Does anyone know how to play a fade out animation when closing the popup?

I am using Leaflet 0.7.3.


texts: Tilemill label and Shield are overlapping, how can this be solved?

texts:
I am trying to use tilemill labels and Shield together for my road designing. But there is a problem with that. Please see the following image for better understanding. The shield is ovarlapping the Road name.

I used the following code for my .mss

#road::Labels {
      text-name: [ROAD_NAME];
      text-face-name: "Candara Bold";
      text-size: 16;
      text-fill: #000000;
      text-placement: line;
      text-halo-fill: #f6f3f3;
      text-halo-radius: 1;
      text-spacing: 400;
      }
    #road::N_Shield {
        [ROAD_TYPE =1]{
        shield-file:  url("C:\wamp\www\shield.png");
        shield-name: [ROAD_NO];
        shield-face-name:"Arial Bold";
        shield-size:16;
        shield-fill: black;
        shield-line-spacing:5;
        shield-min-distance:1;
        shield-allow-overlap:true;
        shield-placement:line;
        shield-dx:150;
        shield-spacing:450;
        }
     }

Anyone can help me what to do.


texts: Why is there a small white space between tiles in Chrome?

texts:
I'm using openlayers to display a google base-map layer and it works fine on Firefox/IE/Safari. I'm having problems with Chrome as the tiles have a small white space between them shown in this image.

I found this discussion about invasive CSS which suggested the problem was  due to

  Your theme is "helpfully" adding margins to img tags. you can override
  this with

.openlayers-map img{ margin:0px }

I'm unclear where I put this line. I read through the openlayers CSS file but I don't have a good idea what is going on, on where I should be editing the img tag.


texts: Get coordinates from a drawing point using OpenLayers 2

texts:
I implemented a drawing feature (OpenLayers2) on an OpenStreetMap.

I searched about how to get coordinates from a drawing point. But I'm still stuck at this point.

I tried this code in order to get the coordinates instantly after drawing a new point but with no success:

 // HOW TO ADD A DRAWING FEATURE
 var DrawMyPoint = new OpenLayers.Control.DrawFeature( MyPointLayer, OpenLayers.Handler.Point)//,
     //{feature:GetPointCoordinates});
     // Point for creating point, polygon creating polygon, and so on.

     map.addControl(DrawMyPoint);

 // EVENT LINKED
 document.getElementById("DrawCheckBox").checked = PushOnButtonDraw();
 document.getElementById("DrawCheckBox").checked =GetPointCoordinates();

 // FUNCTION TO GET COORDINATES OF THE POINT
 function GetPointCoordinates(){
    var coordinates = MyPointLayer.geometry.getBounds();
    console.log(coordinates);
    alert(coordinates);
 }

Here is a link where I saw the information about getBounds or getCoordinates :
Openlayers get coordinates from polygon

I am not sure if I may link two functions with the same event, or is there another solution about it?

With getting the coordinates of each newest drawing point (which is my initial question), I would center the map to this coordinates too.

Do you have any idea?


texts: Add-In button Click open ArcCatalog window with list of .sde files

texts:
I want code to open ArcCatalog window with list of .sde files.The below code always opens an ArcCatalog window with home location. I have to go to database connections from the window which finally gives me list of .sde files. I want code to show list of .sde files directly on click on Add-In button instead of going to multiple selection options in the window and finally listing the .sde files.

                IGxDialog gxd = new GxDialogClass();
                gxd.AllowMultiSelect = false;
                gxd.ButtonCaption = "Add";
                gxd.RememberLocation = false;

                IGxObjectFilter gxObjectFilter = new GxFilterDefaultDatabaseWorkspaces();
                gxd.ObjectFilter = gxObjectFilter;

                IEnumGxObject gxEnumObj;
                gxd.DoModalOpen(ArcMap.Application.hWnd, out gxEnumObj);


texts: Using ogr2ogr to download GeoJSON files from ArcGIS MapServer

texts:
I am trying to download as a GeoJSON file a line feature from a ArcGIS MapServer REST Service run by the National Park Service (link is here) so that I can open the data in QGIS.  Using the guide found here, I typed the following command into ogr2ogr:

 ogr2ogr -f GeoJSON Route66.json "http://imgis.nps.gov/arcgis/rest/services/NHT/NHT/MapServer/9/query?where=objectid+%3D+objectid&outfields=*&f=json" OGRGeoJSON -gt 1000

The command downloaded and created the file just fine (though it did complain about not knowing which CRS to project it to), but when I went to open the file in QGIS it displayed with a number of weird artifacts along the feature.

A picture of what I am seeing can be seen:

Does anyone have any idea how I might go about fixing this situation, if there is any fix to it?  As a side note, I have opened the REST Service data straight from their server on ArcGIS Explorer Desktop, and don't see a problem there.

I'm running QGIS 2.2 from OpenGeo4W on Windows 7 Professional 64-bit.  Thanks for the help in advance!


texts: ArcGIS Server js API - buffer query against polygons and not points

texts:
I'm using the below sample and I want to change the queried results from points to polygons

http://help.arcgis.com/en/webapi/jav...ry_buffer.html

Probably very simple but can't get my head around where to change it.

I tried replacing the queryTask to http://sampleserver1.arcgisonline.co...SA/MapServer/1 and changing the symbol style to reflect this but it didn't seem to make much difference. I'm sure that somewhere and I can't see where, I need to make reference to the fact that we are now querying polygon and not point data.


texts: How to export QGIS map into an interactive map?

texts:
I am using the QGIS 2.18.14 and tried 2 ways to export the QGIS map into an interactive map but failed.

Export the map to a KML format using the MMQGIS Plugin. I thought I can then import this into the Google Map. As the file is too big to directly import to the Google Map. I uploaded to the Google Drive. However, I still could not import the file and received this message An error occurred. You may continue to use the application but any change you make may be lost.
Export the map using the QGIS2Web Plugin. However, even I do not move any of the sub-folder, when I open the index.html, the file is empty.

Any advice how to deal with the above situations?

Or any other way to export the QGIS map into an interactive map?


texts: Google Earth Engine - count non-null values of multiple images and return them as numbers

texts:
I am trying to get quality information for clipped images in an ImageCollection. My approach - so far - was to mask each image and then count the amount of remaining pixels by using "ee.Reducer.count()" and add the result as metadata to each Image.

The problem is that the new metadata "pixel_amount" is cast as a Dictionary and I can't get its numerical values to sort by them or print them as a list.

https://code.earthengine.google.com/42e32bd20eddb3a5a9fb0852306c967c

/////////// GEOMETRY //////////////////////////////////////////////
var potsdam = ee.Geometry.Polygon(
        [[[12.884230720121536, 52.48090537979447],
          [12.884230720121536, 52.2897995867844],
          [13.234419929105911, 52.2897995867844],
          [13.234419929105911, 52.48090537979447]]], null, false);

/////////// CLOUDMASK /////////////////////////////////////////////

function maskclouds(image){
var qa = image.select('BQA');
  var mask_1 =       qa.bitwiseAnd(1 << 4)
                 .or(qa.bitwiseAnd(1 << 6))
                 .or(qa.bitwiseAnd(1 << 8))
                 .or(qa.bitwiseAnd(1 << 10))
                 .or(qa.bitwiseAnd(1 << 12));
  var mask_2 = mask_1.eq(0);
  var mask = ee.Image(mask_2);
  return image.updateMask(mask)
              .select("B.*")
              .copyProperties(image, ["system:time_start"]);
}

/////////// SEARCH DATE /////////////////////////////////////////////

var start = '2016-04-01'
var end = '2016-10-31'
var diff = ee.Date(start).advance(28, 'day').millis().subtract(ee.Date(start).millis());
var list = ee.List.sequence(ee.Date(start).millis(), ee.Date(end).millis(), diff);

/////////// GET IMAGES /////////////////////////////////////////////

function get_images (date) {
  var date2 = ee.Date(date)
var image = ee.ImageCollection('LANDSAT/LC08/C01/T1')
    .filterBounds(potsdam)
    .filterDate(date2, date2.advance(28, 'day'))
    .sort('CLOUD_COVER')
    //.filter(ee.Filter.eq('TARGET_WRS_PATH', 193))
    .map(maskclouds)
    .first()

        return image
}

var LS8 = list.map(get_images)
var LS8_IC = ee.ImageCollection(LS8).limit(10)

/////////// CLIP IMAGES AND COUNT PIXEL /////////////////////////////////////////////

function clip_count_f (image) {
  var clipped = image.clip(potsdam)

    var unmasked_pixels =  clipped.select('B4').reduceRegion({
      reducer: ee.Reducer.count(),
      geometry: potsdam,
      scale: 30,
      maxPixels: 1e9
      });

        var image1 = clipped.set('pixel_amount', unmasked_pixels)

  return image1
}

/////////// CAST TO LIST AND VISUALIZE /////////////////////////////////////////////

var clipped = ee.ImageCollection(LS8_IC.map(clip_count_f)) //.sort('pixel_amount')
var clipped_list = clipped.toList(LS8_IC.size())
print(clipped_list)

Map.setCenter(13.0609, 52.3873, 11)

Map.addLayer(ee.Image(clipped_list.get(6)), {bands:['B5', 'B2', 'B3'], min:8000 ,max:18000}, 'LS8_clipped')
print(clipped)


texts: Custom CRS: different results in PROJ6 and QGIS

texts:
I have created a custom CRS with PROJ6:

echo 115.90428843055600 -1.92307681111111 126.014 | cct +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=tmerc +lat_0=-1.90355333888889 +lon_0=115.900912891667 +y_0=538.374 +x_0=5200.27 +ellps=GRS80 +step +proj=affine +xoff=-5225.30897756133 +yoff=830.070127940344 +step +proj=affine +s21=-0.0072010369955164367159 +s12=0.0072010369955164367159 +step +proj=affine +xoff=5225.30897756133 +yoff=-830.070127940344 +step +proj=affine +xoff=-9.85881089465147 +yoff=-0.172038726321603 +zoff=-18.541

with the results:
 5560.2727     -1623.1404      107.4730           inf

Which is correct.

But when I pass to QGIS (version 3.12.3 with PROJ 6.3.2) the very same proj string I cannot test the values:

If I erase the first +proj=pipeline we obtain incoherent results:

So, how is it possible to obtain the same results in PROJ6 and QGIS (and then used this custom CRS) and why is it not possible, in QGIS, to start a proj string with +roj=pipeline ?

Thanks for your help
Cyril


texts: Color ramp doesn't get applied to categories

texts:
I am following the qgis official training guide, and I'm having this problem. Maybe it's because guide was done for QGIS v. 3.4 while I'm with 3.12, on Debian 10.

I have this vector layer with different classes, and I want different colors for each category in the column "landuse".

To do this, I go into layer properties > Symbology , I use Categorized for Value: landuse and choose a Green color ramp. I classify and what happens is that all categories have the same color, they are all dark green.

EDIT: doing the same exact process with different vector layers does the job just fine

Any idea on what may I be doing wrong here?

This is what I'm trying to achieve:

This is what I'm having

EDIT: this is likely a bug in QGIS 3.12, removing the layer and re-adding it seems to do the trick. It has happen also to some friends on Windows so I don't think is OS related or anything


texts: How to add Oracle database connection to AnalyzeforMSD script?

texts:
I use this script right now:

import os, arcpy

directory = r"C:\Project"
 for root, dirs, files in os.walk(directory):
  for myFile in files:
    fileExt = os.path.splitext(myFile)[1]
    if (fileExt == ".mxd"):
        fullPath = os.path.join(root, myFile)
        print myFile

        myMap = arcpy.mapping.MapDocument(fullPath)
        analysis = arcpy.mapping.AnalyzeForMSD(myMap)

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

It runs on many project mxds. And every time it analyzes a mxd file I`m asked the connection dialogue and have to enter the connection details.

Is there a way to open a database connection that is valid for all mxd?


texts: Manually setting basemap using ArcGIS API for JavaScript?

texts:
Looked over the example can NOT figure out how to set basemap MANUALLY. I don't want a dijit widget, or any other libraries or anything like that. Just want to manually set a basemap to any of the already available types like Topographic, Satellites, Streets, etc. My script code so far is like this.

The part I can't figure out is marked with question marks.

    require([
        "esri/basemaps",
        "esri/map",
        "dojo/domReady!"
    ], function (esriBasemaps, Map) {

    /* ------------------------------------- */
    /* Basemap add one of the existing maps. */
    /* ------------------------------------- */
    esriBasemaps.myBasemap = {
        baseMapLayers ???
    };

    var map = new Map("map", {
        basemap: "myBasemap",
        center: [-118, 34.5],
        zoom: 8
    });

});


texts: Data matching between 2 layers which do not overlap

texts:
I would like to know how a vector layer can takes attributes of a field, of another vector layer which is not overlapping.

For example, there are 2 vector layers, test_1 (2 records A,B) and test_2 (3 records 1,2,3), and i should have as output, the layer test_1 with the same route, but composed/separated by 3 records as test_2, with the additional field ID (1,2,3) moved to the attribute table of test_1.

This geometry process should be as precise as possible.
Any suggestions, QGIS geoprocessing tools?

enter image description here

Thank You for the reply! Using "Join Attributes by Nearest" it doesn't work because the records continue to begins and ends in the same points of Test_1. I should have Test_1 divided in 3 records, for example at the transects of Test_2; with 1,2,3 in the attribute table. Look at the pic.

I'm thinking to use this flow ("Transect" of Test_2 -> "Split by lines" of Test_1 from transect -> "Join Attributes by Nearest" -> "Group geometries" of Test_1 from ID), but this is not accurate when a meet an intersection of lines because the transects create much confusion. Thank you!


texts:
I have several objects showing but they are ordered in my config.xml file with the highest priority on top. Still when I pick a point object that is near a line object it always gives the popup for the line.
I am learning that it is not the order in the mxd that has anything to do with the order in the web display or the published service,
The order in the config.xml does determine the order of the layers (not the sublayers).
1. What determines the order for object selection/mouse pick when there are more than one poup configured?
2. Has anyone successfully manipulated the order of the layers in the service?


texts: FME error when performing a simple process of adding columns to features

texts:
I've created a workbench where I simply need to add four columns into multiple features. I've tried a few different ways. The data has metadata in it so I just want to truncate the table or drop and create so it retains the metadata.

I used a dynamic writer, the error is random.
I tried a more manual way by just exporting "Create if needed" into an empty GDB and importing Metadata afterwards. Still no luck. Weird errors.
I tried an even more manual way by pushing one feature at a time through and having a non-dynamic output. This was random, sometimes it would work, but sometimes it wouldn't. I finally got this to work and pushed through 20 features before it decided to stop working and error "Cannot write a non-point feature to a point feature class" but the output is set to polygon! I've made several new writers and readers and workbenches, and used different features and new GDB's and I can not get this to work.

Please help!


texts: Proper updating of new memory layer geometry?

texts:
Current workflow for a toolbar button:
 -create memory layer
 -save memory layer to shapefile

Then the user uses the built-in QGIS tools to add geometry and finally saves changes via the auto prompt (when done drawing shapes).
When this shapefile is later opened in another QGIS session, no geometry features are present.

I'm inclined to feel the feature addition must occur prior to the memory-to-shapefile save by way of signal slots?
If this is the case, could someone share a simple link or how-to for such?


texts: Force a new AddDataDialog to open at a defined directory path

texts: Is there a way to force the AddDataDialog to open to a specific path?  I created a tool that when clicked on a specific feature will open the Add Data dialog, but I need it to open starting at a specified location.  AddDataDialog doesn't have a property analogous to OpenFileDialog.InitialDirectory.


texts:
I know QGIS comes with its field calculator, which is not bad. It also comes with ftools which permits several operations on geometry and attributes. But we often need more control and I haven't seen any convenient tool in QGIS to perform such operations (where you can select fields you need for example or perform arithmetical operations on them while doing a spatial query).

Basically, I'd like to know if a SQL console (a kind of Mapinfo mapbasic window) is on the roadmap of QGIS and if not, if it would be hard to code it as plugin. If it exists yet, please let me know, i've missed it! ;-)


texts: Creating maps with tile URL using GeoServer?

texts:
I am new to GeoServer, and I plan to create a proprietary map with tile URL.

Will it possible to get the URL like this (http://c.tile.openstreetmap.org/${z}/${x}/${y}.png)
 from GeoServer?

How do I create the maps (If possible)?


texts:
I'm currently having some problem while trying to create a GTiFF file.

require(raster)
require(sp)
require(rgdal)
require(gdalUtils)
library(rgeos)

map <- readOGR ('MyShapefile.shp',layer='MyShapefile')
grid <- paste(sortie,"grid.tif",sep="")
gdal_grid(src_datasource='MyShapefile.shp',dst_filename=grid,zfield='In_ATMO',a="invdist:power=2.0:smoothing=1.0",of="GTiff",ot="Float32",l=mod,output_Raster=T)
sauv<-paste(sortie,"test_In_ATMO.tif",sep="")
gdalwarp(srcfile=grid,dstfile=sauv,tr=c(7.0,7.0),of="GTiff",ot="Float32",overwrite=T)

This is the code sample I use to create my GTiFF file, but it really looks like it isn't in an high definition.

My csv data looks like this

long;lat;In_ATMO
55557754;14454489;4
......

And with this data I create a ShapeFile that I use to create a GTiFF file. if asked, I can share those files to help answering my question..

Is there anyone who knows how can I get a better definition for my Tiff ?


texts:
I am trying to scrape all school boundaries from this website.  I want to download them in any geo-format, such as .geojson, .kml, or .shp

I figured out that the school boundaries are rendered as <path> attributes.  Looking in the source code didn't help.

Does anyone have experience scraping ArcGIS websites, that can give me some hints on how to scrape this ste?


texts: Opening file geodatabase using QGIS for Mac?

texts:
I'm trying to open a GDB on QGIS for Mac but having no luck.

Details:

I'm currently running OS 10.9.3 on my laptop.
I downloaded and installed QGIS 2.6.1-2 from KyngChaos
I installed the GDAL complete package (including NumPy) and the Matplotlib, both from KyngChaos
I attempted to follow Ragi Yaser Burhum's udpated Homebrew advice from this page. I was successful up through step 3; however, step 4 ("Install gdal using brew...") was not successful for me—perhaps because I had already installed it via the package from KyngChaos?

What should I do from here?


texts: How to Select a Polygon beneath another Polygon in QGIS 2.14

texts:
We have a series of layers that contain polygons with additional polygons 'underneath' each polygon.

How do we select the polygon beneath the existing feature without also selecting the feature above it?


texts: QGIS 3.4.3 touch zoom and pan missing from Map Navigation Toolbar is this a permanent loss or will there be a fix?


texts: Is atmospheric correction relevant for Landsat Level2 before classification?

texts:
I am quite new in the field of remote sensing.
I would like to do a separate classification of two Landsat 8 Level2 images and then do a LULCC change detection.
For the pre-processing steps, it is recommended to correct radiometric/atmospheric.
I plan to convert from DN to ground reflectance (or radiance) using atmospheric models (e.g. the i.toar tool from Grass GIS). I have doubts about the relevance of this correction as I am using a Landsat
8 Level 2 image.
Is it mandatory to do this conversion for a L2 image before starting the classification?


texts: How to add and remove the pop-up text based on the condition by using Leaflet

texts:
I want to show and hide the popup text based on the condition given below.
.component.ts
 var insidecirclespace1 = (lat, long) => {
    for (var i = 0; i < this.sensors.length; i++) {

        if (this.getDistanceFromLatLonInKm(this.sensors[i].latitude, this.sensors[i].longitude, lat, long) <= 5) {
            circlemark[i].setStyle({
                color: 'red',
                weight: 1,
                opacity: 7.0
            });

            if (existfalg == 0) {
                this.wifiDrones.push({
                    "frequency": "",
                });
            }
        }
    }

    if (latlngs[latlngidx1] != undefined && latlngs02[latlngidx2] != undefined && latlngs03[latlngidx3] != undefined && latlngs04[latlngidx4] != undefined) {
        outsidespace(latlngs[latlngidx1].lat, latlngs[latlngidx1].lon, latlngs02[latlngidx2].lat, latlngs02[latlngidx2].lon, latlngs03[latlngidx3].lat, latlngs03[latlngidx3].lon, latlngs04[latlngidx4].lat, latlngs04[latlngidx4].lon);
    }
}

Now I want to show the pop-up text when the condition is insidecirclespace1and remove the popup in outsidespace.
I am new to this leaflet map, can anyone help me regarding this.


texts: Constructing distribution graphs of attribute tables using ArcGIS field calculator?

texts:
Since I'm new to python I have an issue with writing the correct code to achieve the following. The picture is an example of an attribute table of which i want to make frequency distribution graphs according to certain ranges of an indicator variable.
For instance a graph should be created in which the range of ORG_STOF (indicator variable) should be: < 2, 2 - 3.5, etc. for which the Count_ORG_STOF needs to be indicated.

According to another question post i found (Using Range in Python expression of ArcGIS Field Calculator?)
I thought that a solution would be something in the line of:
Additional fields should be created for each range (in most cases consisting of 3 or 5)

    def "Range_one"(value):
    if value > 2 and value < 3.5:
     return "Value of Count_ORG"
     else:
     return "0"

Then in the calculation area:

    Range_one(!ORG_STOF!)

However, it is a bit different in my case, since in the solution script I had in mind I need to refer to two fields.


texts: How to only show labels for an arbitrary selection of items?

texts:
I'm curious how others go about solving this problem: You've created a map for something with a large number of features that are labeled. The customer/client asks that you only show labels for X, Y and Z, based on some seemingly arbitrary decision (eg. what they deem important features). How would you go about doing this?

Some ideas:

Create a new string column for this special label and only fill in a value for the features they want to see (could result in duplicate information)
Create a new boolean column and flag the features they want to see with true, then use the conditional labeling in QGIS 1.8 to only display the label when the boolean is true


texts: Adding ArcGIS topology rule that allows Attribute based exceptions?

texts:
I want to run a topology - one is geometry based where I am looking for overlaps as error, but in my database I want to keep some error as exceptions (such as road over waterbodies).
One option is run entire topology and look for exceptions and mark it as exception.

But I have thousands of such exceptions, is there any way to mark them as exceptions based on their attributes?


texts: Change numeric field values into negative

texts:
I would like to change one field values into negative (i.e. original column, row 1 = 5, new column row 1 = -5 in a feature class list.

I got until the addfield and calculate field functions, but it is not clear how to do the negative transformation.

ARCGIS help page has some instructions about negative expressions

-x, negative expression of x,  x = 5 --> -x, -5

How to use it with arcpy is not very clear. It seems that an expression should be added either at the  '!original_value!' or as a code block after the expression type (python 9.3).

for fc in arcpy.ListFeatureClasses():
    #output = path.join(outWS, fc)
    #"GEODESIC" calculates distances in meters!
    #arcpy.analysis.Near(fc, idFeat, "", "NO_LOCATION", "NO_ANGLE","GEODESIC")
    arcpy.AddField_management(fc, 'location_value', "FLOAT")
    arcpy.CalculateField_management(fc, 'location_value', '!original_value!', "PYTHON_9.3", "")


texts: custom tooltip generates "mapView should be present" error

texts:
I am trying to apply the custom tooltip example from http://bl.ocks.org/xavijam/6efe7494eee798a0dbf5 to my map at https://afterfibre.nsrc.org/af_dev.html  The code seems relatively straightforward:

      var aftTooltip = new cdb.geo.ui.Tooltip({
          layer: aftLayer,
          template: '

{{operator}}

',
          width: 200,
          position: 'bottom|right'
      });
      $('#cartodb-map').append(aftTooltip.render().el);

And yet the tooltips do not render and the console produces an error that says:  "Uncaught Error: mapView should be present".  I've googled for this error but not much comes up.  Full code at https://github.com/stevesong/nsrc-afterfibre/blob/master/af_dev.html


texts: How to cut / paste only certain features from a vector layer onto a new vector layer in QGIS?

texts:
I have a large vector layer in QGIS (180,000 features) and I would like to select only certain features (based on their area) then cut and paste them onto a new vector layer. Essentially I want to separate the largest features from the rest of the vector layer so I can work with them independently.

How do I go about doing this?


texts: Detecting line intersections between MultiPolygon and LineString in Turf.js

texts:
I am trying to write a section of JavaScript code that allows a user to analyse a route they have plotted on a map. After the user has finished plotting their route, I want to determine if the route they have plotted crosses over any areas of interest.
To do this I am using the lineIntersect method in Turf.js which should allow me to pass a LineString and MultiPolygon.
I have created a Turf.js MultiPolygon:
{
  "type": "Feature",
  "properties": {},
  "geometry": {
    "type": "MultiPolygon",
    "coordinates": [
      [
        [
          51.50836874890726,
          -0.091989636996398
        ],
        [
          51.50797826040486,
          -0.088994163337176
        ],
        [
          51.507701047593805,
          -0.085230241948645
        ],

(the full object had too many characters to post here so I pasted the full output to CodePen)
and converted the user's route to a Turf.js LineString:
{
  "type": "Feature",
  "properties": {},
  "geometry": {
    "type": "LineString",
    "coordinates": [
      [
        51.50350693825546,
        -0.027165818854444357
      ],
      [
        51.50115610069437,
        -0.041240038776924066
      ]
    ]
  }
}

However when I pass these two objects to the lineIntersect method:
turf.lineIntersect(myLineString, myMultiPolygon);

I am getting the following error:
Uncaught Error: coordinates must contain numbers
    com turf/turf/turf.min.js:1


texts: Fastest way to post layer to GeoDjango map

texts:
I am porting a webmap built with HTML/CSS/Leaflet to a GeoDjango/Leaflet framework, and want to know the fastest way to display map layers. For the HTML/CSS/Leaflet map, my layers are hosted as flat GeoJSON files and called using L.GeoJSON(). The layers display very fast.

For the GeoDjango framework, however, the layers display much more slowly. I use almost identical code, but am guessing the geom data has to be serialized to GeoJSON, which is responsible for the slow display. Is there any way to speed up this process? What would be the fastest way to get layers to display using a GeoDjango framework?

GeoDjango code to display layer:

From urls.py:

from django.conf.urls import include, url
from models import mbls
from djgeojson.views import GeoJSONLayerView

    urlpatterns = [
        url(r'^mbls/', GeoJSONLayerView.as_view(model=mbls), name='mbls'),
    ]

From index.html:

var mblInt = new L.GeoJSON(null, {
                style: { color: '#000000',
                       }
                });

    $.getJSON('{% url "mbls" %}', function (data) {
                mblInt.addData(data);
            });

... other code ...

var map = L.map("map", {
            zoom: 11,
            center: [35.435, -106.54],
            minZoom: 6,
            maxZoom: 20,
            layers: [basemap, mblInt]
        });


texts: QGIS polygonize error 4 "file does not exist in the file system"?

texts:
When I try to convert a raster DEM (in .asc or geotiff format) to polygons I get the error below:

  ERROR 4: `Y:/GIS/oversvmmelser/tmp/out.tif' does not exist in the file system, and is not recognised as a supported dataset name.

Any idea what the problem is?


texts: MapServer pass parameters for a filter through CGI

texts:
Using MapServer and generating a static map.  I am trying to pass a key of a record through CGI and have MapServer draw it.
From what I understand if I have a name=example in my CGI call to my MapServer then inside the MapServer config file I can reference the parameter wit '%name%'. But it does not seem to work.

Here is my layer configruation:

LAYER
    NAME "devlocx"
  GROUP "alarmloc"
    STATUS DEFAULT
    TYPE POINT
    CONNECTIONTYPE postgis
    CONNECTION "host=localhost dbname=osm user=osm password=osm port=5432"
    #DATA "geometry from ptx using srid=3857 using unique gid"
    DATA "geometry from (select gid,name, (st_transform(geometry,3857)) as geometry from ptx  where name='%name%') as foo using unique gid using srid=3857"
    ###########################
    # Define the cluster object
    ###########################
    CLASS
      STYLE
        SIZE 30
        SYMBOL "cir"
        COLOR 255 0 0
      END
    END
END
SYMBOL
 name "cir"
 type ellipse
 points 1 2 end
 filled true
END


texts: Setting different colors based on value when labeling using QGIS

texts:
I want to label my values in a specific color. They are all in the same shapefile which makes it complicated.
For example:

if value is 10, it should appear in blue
if value is 15, it should appear in green
if value is 25, it should appear in red

How can I do that?


texts: Setting FontFamily in OpenLayers Text Label

texts:
Changing the font-Family of a style objects' text label does not have an impact - it always looks the same, no matter what font-family I set. Do I have to provide the respective fonts to openlayers somehow or do I have to do something else?

            var fontFamily = "Arial"; // "Verdana" // "Courier New";
            style.setText(new ol.style.Text({
                text: 'my-label',
                font: '12px bold '+fontFamily,
                overflow: true,
            }));

Funnily, fontFamily = "Helvetica, sans-serif"; works, while fontFamily = "Helvetica"; and fontFamily = "sans-serif"; do not work (i.e. is displayed with the default serify font). Any explanation for that?


texts: How to serve local maps for an offline HTML5 webmapping app?

texts: How could one serve maps (tiles in wherever format possible) within a Android device to be consumed by an offline webmapping app also accessed within the device's browser?


texts: Creating descriptor out of other Descriptors Using ArcGIS Network Analyst

texts:
I created a Network Analyst route, and I just added a few descriptors in it.
I want to know if there is a way of create a descriptor based on other descriptor.

The equation is :

Descriptor 4 = Descriptor 1 / (Descriptor 2 + Descriptor 3)

Is that possible to do on ArcGIS Desktop?


texts: ogr2ogr reprojection works on command line but shifts output when same command run from inside python?

texts:
Reprojection with ogr2ogr works fine on command line, but shifts the polygon to northeast when same command is run from inside python using os.system(...)

Any ideas why this could be?

I hard-coded the full ogr2ogr path to avoid hitting different installations from each environment.

PRJ files for both versions are identical so the difference must be in the "shp" file.

e.g. Here is the command:

/opt/local/bin/ogr2ogr -s_srs epsg:28992 -t_srs epsg:4326 -a_srs epsg:4326 /mnt/output_4326.shp /mnt/input_28992.shp


texts:
I would like to use QGIS to convert a CSV file from UTM to lat/long.

I imported the CSV file with the "Add Delimited Text Layer" plugin.  This works fine, the data looks like it is projected properly on my map.

However, when I try Layer> Save as... and set the Selected CRS to NAD83, it does not convert my UTMs to lat/long.  When I go Open Attribute Table, the x and y coordinates are still in meters, not degrees.

Any thoughts why this isn't working for me?


texts: Export data with data frame coordinate system using ArcPy

texts:
I want to export data in current data frame and use the same coordiante system as data frame using ArcPy as I have multiple files. Any suggestion on how to "Export with data frame's coordinate system"? Following is what I have:

import arcpy
from arcpy.sa import *

mxd = arcpy.mapping.MapDocument(r'C:/test.mxd')
df = arcpy.mapping.ListDataFrames(mxd, 'Layers')[0]
for layers in xrange(1,12,1):
    inLU = r'C:/rl_' + str(layers) + '.lyr'
    lyrFile = arcpy.mapping.Layer(inLU)
    arcpy.mapping.AddLayer(df, lyrFile)
    arcpy.CopyFeatures_management(r'rl_' + str(layers), r'C:/shape_' + str(layers))


texts: Fitting an ellipse within an irregular polygon shapefile

texts: Can we fit ellipse in irregular polygon shapefile in QGIS like Find maximum radius of circle that will fit within an irregular polygon?


texts:
In ArcGIS ModelBuilder, is it possible to create a variable (that can be used as a precondition for following tools) that tells if an iteration has been finished?

Example:
An iterator buffers all feature-classes in a file-geodatabase and creates an output for each of that feature-classes. After all feature-classes have been buffered, I want to run a new tool.

How do I set the completed iteration as a precondition for that tool?

I'm not fluent in python.


texts:
I have a basemap w/ GPS points on it, and polygons of certain areas.  The polygons fall well within the borders of the basemap.  I am trying to generate a Euclidean Distance map of distance to these polygons, with the ENTIRE basemap as the Processing Environment.  The problem I keep running into is that the Euclidean Distance generates a map that only covers the absolute (x,y) extent of the polygon files.  I have tried setting the processing extent as the entire basemap, converting the polygon file to a raster with the basemap as the extent.....to no avail.

How do I generate a Euclidean distance map that goes beyond the input feature extent?


texts: I want to know slope value between two point with lat and long coordinates. I can also get DEM values at lat long location for both points. What should I do next in order to calculate slope?


texts: Extracting a Feature Class' Name, and Adding It to a Field

texts:
I'm building a bike network dataset based on an existing street network. Each bike lane/route/path has a unique ID number (assigned by our engineering deparment) that must be included in a separate ID field within each feature class. We recently changed the base for our street network, necessitating a redraw of about 4000 segments.

To save time manually adding fields, I built a model to add the new ID field "BPID" (which must be saved as a short integer), a "BPID_Temp" field, used to convert the string values to numbers by trimming the leading character (expression is =Right( [BPID_TEMP], Len( [BPID_TEMP] )-1 )  and then copying it to the permanent BPID field. Lastly, I delete the temporary ID field.

Here is an image of the model below:

At the first Calculate Field step, the value from Parse Path is not being transferred into the BPID_TEMP field. The expression for the step is "Value", in PYTHON. Parse Path is set to NAME, and the model runs through without throwing any errors.

I've gone through Using Feature Class "Name" from Iterate Feature Classes to Calculate Field with ModelBuilder?. However, adding wildcards to the "Value" expression returns an invalid expression error. If I link the "Value" output from Parse Path directly into the Calculate Field (rather than as a precondition), it tells me that I am trying to calculate a number that is too large, or mixing incompatible data types.

Am i missing something critical here. I can run the steps manually no problem, but as soon as I try to automate with an iterator, it all falls apart. I'm not great with Python yet, but I feel like it isn't even necessary here.


texts: How to use union in area queries (QL/XML) with Overpass API/OpenStreetMaps?

texts:
I am using area queries to find nodes with particular amenities Eg:-

[out:json];
area[name = "New York"];
node(area)[amenity=bar];
out;

However the issue with this is that if I query the following:

[out:json];
area[name = "New York"];
(
  node(area)[amenity=bar];
  node(area)[amenity=restaurant];
);
out;

The output includes only the first amenity, i.e. 'bar' even though there's a union. If I reverse the order to put 'restaurant' before 'bar' the output includes only restaurants. I also tried the same using XML syntax but in vain. Any ideas?