texts: Has Action behaviour changed between 1.7.4 and 1.8.0?

texts:
I have an Action that loads raster layers (.jp2) into the current QGIS window.  This action has been working great in version 1.7.4.  I have now upgraded to 1.8.0 for many reasons, but the same Action is not working.  Here is the Action I have been using: qgis.utils.iface.addRasterLayer("E:\\Projects\\CityTBay\\Imagery\\2008\\%image")

The Action will load an ortho image from the path location specified.  The image name is saved in the filed image.  When I use this Action in 1.8.0 the following error occurs:
E:\Projects\CityTBay\Imagery\2008\%image is not a valid or recognized raster data source.

The Action is using the Python option.

Thanks...


texts: Coordinates must be finite numbers error on map.getBounds()

texts:
I get the following error when trying to get a getBounds for my map.
*proj4-src.js:1849 Uncaught TypeError: coordinates must be finite numbers
    at checkCoord (proj4-src.js:1849)
    at checkSanity (proj4-src.js:1841)
    at transform (proj4-src.js:1865)
    at transformer (proj4-src.js:1944)
    at Object.forward (proj4-src.js:2010)
    at NewClass.project (proj4leaflet.js:40)
    at vers copy.html:80*

I'm using a 28992 projection.
  var RD = new L.Proj.CRS(
        'EPSG:28992',
        '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812 +no_defs', {
        origin: [-285401.92, 22598.08],
        resolutions: [3440.640, 1720.320, 860.160, 430.080, 215.040, 107.520, 53.760, 26.880, 13.440, 6.720, 3.360, 1.680, 0.840, 0.420, 0.210, 0.105],
        bounds: L.bounds([-285401.92, 22598.08], [595401.920, 903401.920])

var wgs84 = map.getBounds().getEast();
console.log(wgs84);
RdNew= RD.projection.project(wgs84);
console.log(RdNew);

The console log for wgs84 gives the correct coordinate.
The console log for RdNew gives the error listed above


texts: Dealing with zeros while iterating a model for percent change

texts:
I'm using ArcGIS to examine some aquatic vegetation changes over time and I have created a model that adds two fields: Area_percentLoss and Biomass_percentLoss. My model is iterating over feature classes in a gdb (annual vegetation polygons). My calculate field for percent of biomass lost is

([BiomassLoss_kg]/[TotalBiomass_kg]) * 100

However, when it comes to a BiomassLoss = 0 I get an error. How could I write a Python script that does something to the effect of:

Biomass_PercentLoss =

IF BiomassLoss = 0
THEN Biomass_PercentLoss = 0
ELSE Biomass_PercentLoss = ([BiomassLoss_kg]/[TotalBiomass_kg]) * 100

It's just dividing something by 0 is throwing an error, I think.


texts: flex 2.3.1 Simpleinfowinwidget deprecation

texts:
I have been using the simpleinfopopup widget and recently read that it is being deprecated with 2.3.1.
1. What are the positives and negatives of the popup you are using in flex?
I also have been trying to implement the link_field without success in this popup.
2. Have you had any success getting a url link in your popups?


texts:
I try to download data from the city of Rome with this command:

ogr2ogr -f GeoJSON EVC_rome.json "http://services2.arcgis.com/NZMqCJwY3kMjFOqf/arcgis/rest/services/ricarica_elettrica/FeatureServer/0/query?where=objectid+%3D+objectid&outfields=*&f=json" OGRGeoJSON -gt 1000

The OSGeo4Q Shell says "unable to open datasource `http://services2.arcgis.com/NZMqCJwY3kMjFOqf/arcgis/rest/services/ricarica_elettrica/FeatureServer/0/query?where=objectid+%3D+objectid&outfields=*&f=json' with the following drivers."

Any ideas?


texts: Saving Options file for ArcGIS Desktop to release hanging licenses?

texts:
I have created an options file for ArcGIS 10.1 so that I will be able to release a hanging ArcGIS license after a set amount of inactivity. I found that for ArcGIS 10.0 the correct directory to store the Options file on is C:\program files\ArcGIS\License10.0\bin.

Unfortunately for ArcGIS 10.1 once I open up the ArcGIS folder in C:\Program Files\ I can not figure out where to save the Options file.


texts:
I'm making an addin for ArcMap with ArcObjects and I want to detect when the user uses the split tool on a line feature. How can I achieve this?

    public StartUp()
    {
        var editor = MapHelpers.GetEditor(); //Returns the application's IEditor3
        var events = (editor as IEditEvents_Event);
        events.OnChangeFeature += (obj) =>
        {
        };
    }

So far all I can do is hook into the OnChangeFeature event and get the resulting IObject. I don't know where to go from here as the documentation is lacking.


texts: Create ArcGIS geodatabase domains from a CSV file using Python

texts:
I want to create ArcGIS geodatabase domains from a CSV file.The problem is that my output says field values in the row are "None". It doesn't recognize the values in my CSV file.
Any ideas ?
Here's the code I am using :
import arcpy
import csv
arcpy.env.workspace = "D:/DATA"
gdb                 = "NEW.sde"
csvfile             = "D:/DATA/Domains.csv"
domains = []
rows = []

with open(csvfile, 'r') as csvfile:
    reader = csv.DictReader(csvfile)
    for row in reader:
        try:
           domains.append(row.get('domain'))
           rows.append(row)
        except KeyError:
           rows = [row]

 domains = list(set(domains))
 print("Identified %s unique domains" % len(domains))
 print("Identified %s codes" % len(rows))

 print("Creating codes:")
 for row in rows:
     print("%s\t%s\t%s" % (row.get('domain'), row.get('code'), row.get('code_description')))
     arcpy.AddCodedValueToDomain_management(gdb, row['domain'], row['code'],  row['code_description'])

print("Done")

And here's my CSV File:
[![enter code here][1]][1]


texts: Trying to run ArcMap says appropriate license not given to program?

texts:
I am trying to run ArcMap 10.4 on Windows 10 with a checked out license, and I get this error:

When I turn off my Wi-Fi and intiate ArcMap, it works, then I can turn the Wi-Fi back on and continue to do what I need to.

Have any of you encountered this?


texts: Locating Coordinate Reference System Selector on QGIS

texts:
I downloaded the latest version of GQIS last Friday(6/3/15). I can't locate the Coordinate Reference System Selector.
I have tried clicking settings, custom CRS, add CRS, but can't find a way to proceed from there. Any help?


texts: Find shortest path between two points without spatial database

texts:
I have an an certain area contains less than 200 roads(polylines), and I need to calculate a shortest path between two specified points inside the area.

I know maybe I can use the pgroute in postgis to do the calculation, however I wonder if there is an alternative to use common database like mysql instead of spatial database since the dataset is not large enough?


texts: Adding an AddIns.Button to an Add-In Menu

texts:
I know the title sounds unwieldy but as I am beginning convert an existing COM interface to Add-Ins, I am having trouble adding a button to a root menu.

I have an AddIn.Button that is used both as a button (in Add-In Controls) and an item in a rootmenu. The button works but when I add the namespace.classname to the .esriaddinx (as a Button refID="...", I get the modifier "static' and "internal" errors from the config.designer.cs. How do I reference this addin.button in the rootmenu so it can build the designer correctly?


texts: PostGIS 2.0 geographic polygon has max longitudinal range

texts:
My problem is that I have a load of geographic point data (srid 4326) in a table.  A user may do a bounding box search - essentially create a polygon as follows:

l_poly geography(Polygon,4326):= 'POLYGON((-120 -89, -120 81, 59.999 81, 59.999 -20, -120 -20))';

... then query using the ST_Intersect predicate

WHERE ST_Intersects(my_tab.point_data, l_poly)

In the above example, I can return a point at lon=0, lat=0.  However, once my eastern-most longitude goes above 60 degrees (currently 59.999), then the longitudinal range exceeds 180 degrees, and my polygon wraps the opposite way around the earth.  In this case I can not return data at 0 longitude, but I would at e.g. 180 longitude.

It seems as though a geographic polygon cannot be defined that exceeds a longitudinal extent of 180 degrees.  Is this true?  Can anyone point me in the direction of a solution?

I have tried casting a geography from a geometry defined with enforced Right Handedness, which I thought might define the interior of the polygon to the right hand side of a clockwise drawn polygon.  But this doesn't work either.


texts: Assign raster values based on polygon intersection

texts:
I have a raster where every cell is a 1. I also have a polygon of buffer areas that fall within the raster. I want to assign all raster cells that intersect the polygon geometrics to NA. But the two files have different extents and origins. I don't see a reason to share code, since I build the raster and polygons based on proprietary data. But I can illustrate it. The black grid represents my current raster, the blue cicle is the polygon:

What function will allow me to reassign the raster value of all cells that intersect the polygon? The raster would become as follows:


texts: ogr2ogr gdb to postgis -- empty geometry after import

texts:
I am using the below command:

ogr2ogr -f "PostgreSQL" -s_srs EPSG:2286 -t_srs EPSG:2286 PG:"host=${PGHOST}
dbname=${PGDATABASE} user=${PGUSER}" import_file.gdb -overwrite -lco
SCHEMA=some_schema -lco OVERWRITE=YES

The import runs without erros and I end up with features having a column

wkb_geometry geometry(Point,2286)

but there is never anything in it.
When I open the same file with FME Safe, I can see the points on the map, so I'm assuming the data has geometry info, but for some reason I'm not importing it?
I've imported gdb dumps in similar manner (from different sources) in the past and never had this issue before..

The output of

 ogrinfo -al -so some_dump.gdb

(ogrinfo -al some_dump.gdb) produces the same exact output

Layer name: SOME_LAYER_FEATURE
Geometry: Point
Feature Count: 344
Extent: (1787446.125000, 496713.522000) - (2007796.551000, 958150.241000)
Layer SRS WKT:
PROJCS["NAD_1983_StatePlane_Washington_South_FIPS_4602_Feet",
    GEOGCS["GCS_North_American_1983",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.0174532925199432955]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["False_Easting",1640416.666666667],
    PARAMETER["False_Northing",0],
    PARAMETER["Central_Meridian",-120.5],
    PARAMETER["Standard_Parallel_1",45.83333333333334],
    PARAMETER["Standard_Parallel_2",47.33333333333334],
    PARAMETER["Latitude_Of_Origin",45.33333333333334],
    UNIT["Foot_US",0.304800609601219241]]
FID Column = OBJECTID
Geometry Column = SHAPE
FID: Real (0.0)
ORIENTATION: Real (0.0)
Z: Real (0.0)

I did just noticed that the line features DO contain geometries.. but the points do not.. even though I can see the points when opening the dump from FME.
I know the dump was generated through FME if that makes any difference..


texts: How to convert date to an arbitrary format in FME Workbench?

texts:
My source dataset is DWG and destination is GDB. The source dataset has Date attributes as Char in DD-MM-YYYY. The corresponding destination attributes have to be of Date datatype, that too in M/D/YYYY format.

I tried DateFormatter transformer, with output date format string as FME Date (%D), I got the error as

The date 'FME Date (05/14/1981)' is not valid. Valid forms for a date type are YYYYMMDD, YYYYMMDDHHMMSS, or HHMMSS

For format string as %D, I got error as:

The date '05/14/1981' is not valid. Valid forms for a date type are YYYYMMDD, YYYYMMDDHHMMSS, or HHMMSS

For format string as %m"/"%e"/"%Y and %m/%e/%Y, again the above error message.

As from the error, I cannot perform the desired formation. Is there a transformer available in FME Workbench to perform desired format conversion? Or should I try something like using SubstringExtractor and then StringConcatenator? I need to read the documentation thoroughly to see if I can get it working.

Thanks


texts: Algorithm to find distance of an uneven polygonal road

texts:
While working for my summer internship I have got a task to implement an algorithm to find the crosssecting points of two roads or more and then place that cross-section point a number and then find the distance of the previous roads to the new point and update it in the table.

I tried Ramer Douglas Peucker Line Simplification Algorithm to convert the road into a straight line and then find the intersection point but the point coming after this is not exactly the cross section point and the distance coming is also not the exact.

Please help.


texts: Modelbuilder iterator only puts out the last input?

texts:
I need to use an iterator in a model process. My problem is that I always only have one output instead of three.
There are three Polygons e to run. The output always only contains the last Polygon one and I have no clue why. I did name the output xy* Value * so that cannot be the problem.

Anyone any ideas?

I'm gonna explain it a bit more detailed.
I've seen the other question with exactly the same title but I couldnt find an answer there.

I need do build a cost distance model which is based on a shp with 3 startpoints as input fpr cost distance tool and 3 shp with each containing an goal point as input for the cost path tool. I want to use an iterator so that it can run all three points with one run. The model should also be with variables so that it can be used for many regions. This is how it looks like with variables:

The names are written in German but I guess you can understand it. It works perfectly exept of the fact that it only puts out the last cost path. So the iterator works and runs each path correctly but only saves the last one.
I did name the output with %Value% (thats what I meant yesterday with Value.. I was just a bit too tired to realise that I didnt wright it correctly..).
So that cannot be the reason why it doesnt save the first two outputs.

I've also tried to run it with given inputs without variables but that didnt make it any better:

I also already tried collect value but that didnt change anything.

I really hope you guys have an idea.

Another problem is that I dont understand how I have to use an iterator with variables. Because the output name is meant to be empty. That means I cannot write %Value%.


texts: Masking Landsat-8 and Modis (Aqua) with same shape file doesn't result in same shape

texts:
I was trying to mask Landsat-8 and Modis (Aqua image). FYI, I have  converted the MODIS image to the same projection as Landsat and same resolution as 30m.  Finally once i am done with masking using the shape file the result comes out as shown below in QGIS. I have two questions regarding the same

1) As you can see width and height of both the output is not same. It differs by 1. I have tried using with different shape file also but again the issue remains same. I always get height 1 less than MODIS. Though i also tried setting the value of nodata to be 0. But no luck.

2) As you can see below both the images are not overlapping perfectly in QGIS screenshot. Moreover, their extent is also different. So, how do i overcome this problem?

Landsat:

Modis

Overlay of Modis on top of Landsat


texts: What are license options to share and distribute ArcGIS ModelBuilder toolkit?

texts:
Do you have expertise or experience to share on what software license can be used to share a model created in ArcGIS ModelBuilder to ensure that it is free to use, or used for non commercial purposes.

I am currently looking at the GPL v3 license, but since the model is based on using ArcGIS software I am not sure if this is appropriate.

UPDATE

I should have added for clarity.
We wish that although the Models are free to use that there is the ability to restrict commercial use (which may then be under a different license)
We wish to be able to ensure that any edits or improvements to the Models are also made available and free to use.


texts:
I have a DEM that was made from UAV/photogrammetry techniques by Maps Made easy.
In QGIS I loaded the DEM and then made a slope raster layer.
My values go from 0 to a huge number like 23456e+07. Does that seem right?
The highest slope on the terrain is not much more than 28%.
I think the layer is spitting out wrong data. Do I need to enter a Z factor?
Also, when I enter my formula in raster calculator ("slope@1" > 28) * 0 + ("slope@1" < 28) * 1 it turns my whole image black.

My end goal is to display a map with the georeferenced jpeg with areas over 28% slope shaded in red.


texts: Harmonizing Data Across Geographic Scopes/Layers

texts:
A common problem in the analysis of spatial data is deciding how to harmonize data from different levels of geographic resolution.  For example, suppose Municipality $X$ crosses a county border, so that a portion of $X$ lies in County $A$, while another lies in County $B$.  If the municipality is my unit of analysis, but detailed socioeconomic data exist only at the county level, how do I map that data to $X$?  The two approaches that I have seen are land area-based proportional allocation and plurality allocation.  I find both of these unsatisfying, so this question is an effort to discover other approaches.

Land Area-Based Proportional Allocation

In this scheme, let $p_A$ be the proportion of County $A$ that is taken up by $X$, and $p_B$ be the parallel proportion of County $B$ taken up by $X$.  Assume that structures and people are uniformly distributed across both counties.  If $n_A$ and $n_B$ are the population counts in County $A$ and County $B$, respectively, the population of $X$ ($n_X$) is given by $p_A n_A + p_B n_B$.

This approach has a nice intuitive feel insofar as it views $X$ as a composite construct with characteristics from both counties.  However, the uniform distribution assumption is a heroic one, and if we want to speak about anything other than simple counts (e.g. the distribution of income) even stronger assumptions are needed.  (One could envision an assumption of randomly distributed incomes based on the distributional information in each county, and then folding the distributions of $A$ and $B$ together.)

Plurality Allocation

In this scheme, let $p_{X \in A}$ and $p_{X \in B}$ be the proportions of $X$ that reside in County $A$ and County $B$, respectively.  Suppose $p_{X \in A}>p_{X \in B}$, which in this case means the majority of $X$ resides in County $A$ (a plurality suffices when more than two counties are involved).  Municipality $X$, then, takes on all the characteristics of County $A$, and County $B$ is forsaken.

This approach avoids the cross-tab pitfalls of the proportional allocation approach.  It does, however, fail a face validity test when $p_{X \in A}$ is low.  For example, if $p_{X \in A} = 0.51$, do we really believe that $X$ is appropriately represented with County $A$'s characteristics?  Furthermore, if I am simulating an activity that relies on this method, how does one true up the simulated aggregates at the state level?

Ideas?

There must be a better way.  Any suggestions?


texts: Line defined by two geographical points crossing circle area with defined radius

texts:
I'm trying to calculate the number of lines crossing a circle area defined by a center and radius.

This algorithm didn't work for me, it returns too many lines.

Is there any other approach which could help my to calculate the lines crossing the circle?


texts:
I have a virtual layer in QGIS 2 with a query on a CSV. The CSV gets updated regulary. I need the max/min attributes in the virtual layer to be updated as well. The virtual layer is used as a title in print composer (displayed as attribute table).

The query looks like SELECT "Events between " as from, min(fromdate) as fromdate, " and " as to, max(todate) as todate from "Events_Happening" WHERE invalid = 0
(Events_Happening is the .csv)

How can I achieve this without recreating the virtual layer by hand every time the csv is updated? (Update on QGIS startup would be enough)


texts: Dynamic loading of KML files into Google Maps

texts:
I'm looking to load KML files that I have created in to Google Maps. I have loaded them to dropbox to test and made public. Can anyone see what I'm doing wrong in the following call:

http://maps.google.com/maps?q=https://www.dropbox.com/s/kgpbcqdnpgqpr2q/Fraud_Google_API_MAP.kml?dl=0"


texts: Where can I find the version of scipy python package used in QGIS Lyon 2.12?

texts:
I am using QGIS Pisa 2.10. I want to use the Voronio class of scipy.spatial in QGIS python console. QGIS Pisa has 0.11 version of scipy package. This package doesn't have Voronoi class. I even tried upgrading it. But since I am in windows environment, the Scipy.exe installer is unable to detect the python inside Pisa (C:\Program Files (x86)\QGIS Pisa\apps\Python27). Instead it is always referring to the Python installed at the user level (C:\Python27_32Bit) from the registry.

If the version of scipy used in QGIS Lyon is greater than 0.11 then, I will uninstall Pisa and install the Lyon.


texts: Passing list containing file paths into function arcpy.Exists to determine if shapefile exists in directory

texts:
I have a dBase table that has a field containing file paths. Within this field there are over 1,000 records. I created a list, using a search cursor, that contains all the file paths. An example of an element in the list is C:\TestFolder\Folder1\A\address_pt.shp

I want to pass this list to the arcpy.Exist function one by one to see if the shapefile exists within the defined workspace.

The code that I currently have only creates a list and prints out the items in the list.

import arcpy
arcpy.env.workspace = r"define directory for arcpy.Exist to look in"

#Set variables:

myLayer = "define the dbf path"
myField = "APFILE_PATH"

myList =[row[0] for row in arcpy.da.SearchCursor(myLayer,myField)]
#print(myList)

for item in myList:
    print (item)
#Insert Arcpy.Exist somewhere in this for loop?

I am new to python and arcpy and I'm not sure how to structure the for loop to take an element out of the list, one by one, and insert it into the arcpy.Exist function to see if the shapefile exits.

Ideally, I want to see if the shapefile exists, based on the file path, and if it exists I want to mark it as a Yes or a No, possibly within a new field within the existing dbf file. Also, if the shapefile exists, I will eventually pass it along to another function that will check if fields exist within the shapefile. However, for now, i'm just working on understanding and coding the first part to see if the shapefile exists.

I'm not sure if this is even possible with the arcpy.Exist function. I'm taking this process step by step.


texts: What is (if there's any) difference between MAUP and Edge Effect

texts:
I heard both the term MAUP and Edge Correction in association with Spatial Statistics, zoning/quadrats, demographics etc.

Is there a difference between the two? What is the difference then?


texts: How to record an unknown geometry in QGIS?

texts:
In the QGIS Python API  it states that "a null geometry will always return true for isEmpty()."

For a geometry, a relational database null refers to an unknown geometry.  For example, if I had a database table of roads, for a particular road I might know the name and even the suburb, but, I don't know the line (representing the path the road takes).  This could happen because one person in my team enters the text attributes and then another digitises the road.  In this example, if I don't know the line for the road path I give a null (meaning unknown).

Now, say I have a database table that has one column for suburb boundaries (polygon), and, another column contains the lines of all the roads in each suburb (multi-line).  Now, when I have a suburb without any roads then I put an empty geometry in the line containing column.  This empty geometry lets me know that it is known that that there are no roads (lines) in a particular suburb (polygon).

The logic that we use in a database like PostGIS is that empty is known and one of many possible values for a null as it is unknown.  But, in QGIS the logic is that null is one of several possible values for empty.

Is there a test in QGIS for an unknown geometry that does not imply that the geometry is also empty?


texts: Point cluster symbology for different categories without overlapping symbols

texts:
I'm trying to visualize crimes based on if the crime happened during the day or night time. I've tried using the point displacement tool but I'm not satisfied with the results because the readability for certain scale levels is too low and the results also overlap:

So instead I followed Anita Graser's guide and tried using the point clustering-tool but I can't think of a way to separate day/night crimes without duplicating the layer and setting different offsets:

Ultimately I want to show the clusters including the exact amount of crimes separated by day and/night without having overlapping feature symbology. So I tried to accomplish this for one category by using the below expression but the results are not correct, the total count is not aggregating correctly:
CASE
     WHEN "day_night" = "day" THEN @cluster_size
END

Using point cluster symbology in one layer, how can I re-write the expression above to correctly count the numbers of day/nigh crimes separately and visualize it without overlaps?


texts: Getting x-coordinate of midpoint from line layer

texts:
I Need to get a midpoint from a line layer.
Property layer -> Style -> Simple symbol -> Generator geometry (line)
I have code #1:
x(line_interpolate_point($geometry, $length / 2)

this code return NULL
I have code #2:
point = $geometry.interpolate($geometry.length() / 2).asPoint()
make_line(make_point(point.x(), point.y()),
make_point( "auxiliary_storage_labeling_positionx", "auxiliary_storage_labeling_positiony"))

this code return NULL
What is my mistake? I am using QGIS 3.8.3
Answer eurojam:
Edit code result:
 make_line(
 make_point(x(transform(line_interpolate_point(transform($geometry,'EPSG: my geographic','EPSG: my project')),'EPSG: my project','EPSG: my geographic')), y(transform(line_interpolate_point(transform($geometry,'EPSG: my geographic','EPSG: my project')),'EPSG: my project','EPSG: my geographic'))),
 make_point( "auxiliary_storage_labeling_positionx",
 "auxiliary_storage_labeling_positiony")
 )


texts:
I am trying to automate the export of NetCDF files into GeoTIFF. However the export could not give me the same result that I do it manually. I am quite new so wasn't sure which part went wrong.
This is the information of export I wanted:
Information from provider
Name    test_compare
Path    D:\WAIG\Progs\qgis_nc\compare_output\test_compare.tif
CRS EPSG:32649 - WGS 84 / UTM zone 49N - Projected
Extent  152123.7194000000017695,85660.8351999999722466 : 847876.2806000001728535,611607.4255000000121072
Unit    meters
Width   70
Height  53
Data type   Float32 - Thirty two bit floating point
GDAL Driver Description GTiff
GDAL Driver Metadata    GeoTIFF
Dataset Description D:\WAIG\Progs\qgis_nc\compare_output\test_compare.tif
Compression
Band 1
STATISTICS_MAXIMUM=100984.875
STATISTICS_MEAN=100244.48617124
STATISTICS_MINIMUM=91913.625
STATISTICS_STDDEV=1502.4594772555
STATISTICS_VALID_PERCENT=100
More information
AREA_OR_POINT=Area
Dimensions  X: 70 Y: 53 Bands: 1
Origin  152124,611607
Pixel Size  9939.322302857144678,-9923.520571698114509

Where the result I get from the script is:
Information from provider
Name    testing_written
Path    D:\WAIG\Progs\qgis_nc\test_output\testing_written.tif
CRS EPSG:32649 - WGS 84 / UTM zone 49N - Projected
Extent  152123.7194000000017695,85660.8351999999722466 : 847876.2806000001728535,611607.4255000000121072
Unit    meters
Width   70
Height  53
Data type   Float32 - Thirty two bit floating point
GDAL Driver Description GTiff
GDAL Driver Metadata    GeoTIFF
Dataset Description D:\WAIG\Progs\qgis_nc\test_output\testing_written.tif
Compression
Band 1
STATISTICS_VALID_PERCENT=0
More information
AREA_OR_POINT=Area
Dimensions  X: 70 Y: 53 Bands: 1
Origin  152124,611607
Pixel Size  9939.322302857144678,-9923.520571698114509

I am not sure how to export the bands together as well, I have read the documentation but at my stage its quite hard to understand... But anyway here's the code:
from qgis.core import *

#use path as input
uri='NETCDF:"D:\WAIG\Progs\ERA5_download\Downloaded_nc_file\sample_input/era5_surface_pressure_2009-12-27_0000.nc":sp'

#create layer to store nc file
rlayer = QgsRasterLayer (uri,'layer20091227_0000')
rlayer.setCrs ( QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.EpsgCrsId))

pipe = QgsRasterPipe()
pipe.set(rlayer.dataProvider().clone())

file_writer = QgsRasterFileWriter('D:\\WAIG\\Progs\\qgis_nc\\test_output/'+'testing_written.tif')
file_writer.writeRaster(pipe, 70, -1, QgsRectangle(152123.7194000000017695,85660.8351999999722466,847876.2806000001728535,611607.4255000000121072),  QgsCoordinateReferenceSystem(32649))


texts: Are Listing data functions of ArcPy correctly documented with respect to their Return Value?

texts:
ArcPy has a number of functions for listing data, and I will use just one of them as an example.  However, what I ask seems to apply to all of them.

In the Syntax section of the help for ListFeatureClasses it says that the Data Type of its Return Value is a String.

However, whenever I use the function what is returned is a Python list, which is what I would expect and desire.

Is there some convention that I am not aware of whereby a Python list can be classified as a string, or should this Return Value be documented as a List instead?


texts: How to generate a collection of mosaic images by mapping over date list in Google Earth Engine

texts:
I want to generate a collection of mosaicked images from Landsat 5 over multiple years. I have checked that the individual images are not null but on returning the mosaic in the code below throws the error:

  List (Error) ImageCollection.fromImages: Attempt to create an ImageCollection with non-image elements.

var p1 = ee.Geometry.Point([103.521,13.028]);
var p2 = ee.Geometry.Point([105.622,13.050]);

var Date_Start = ee.Date('2000-05-01');
var Date_End = ee.Date('2007-12-01');
var Date_window = ee.Number(30);

// Create list of dates for time series
var n_months = Date_End.difference(Date_Start,'month').round();
var dates = ee.List.sequence(0,n_months,1);
var make_datelist = function(n) {
  return Date_Start.advance(n,'month');
};
dates = dates.map(make_datelist);
var fnc = function(d1) {
  var start = ee.Date(d1);
  var end = ee.Date(d1).advance(14,'month');
  var date_range = ee.DateRange(start,end);
  var S1 = ee.ImageCollection('LANDSAT/LT5_L1T_TOA')
    .filterDate(date_range)
    .sort('CLOUD_COVER')
    .filterBounds(p1).first()

  var S2 = ee.ImageCollection('LANDSAT/LT5_L1T_TOA')
    .filterDate(date_range)
    .sort('CLOUD_COVER')
    .filterBounds(p2).first()

  var mosaic = ee.ImageCollection([ee.Image(S1), ee.Image(S2)]).mosaic();

  return mosaic
};

var list_of_images = dates.map(fnc);
print('list_of_images', list_of_images);
var mt = ee.ImageCollection(list_of_images);
print(mt);
Map.addLayer(mt, {}, 'mt');

ERROR:

List (Error)
    ImageCollection.fromImages: Attempt to create an ImageCollection with non-image elements.


texts: Searching an R alternative to ArcGIS Combine tool (Spatial Analyst)

texts:
I'm looking for a methode in R to do the same thing as the ArcGIS Combine tool (Spatial Analyst).

The tool is explained here: https://pro.arcgis.com/fr/pro-app/tool-reference/spatial-analyst/combine.htm

I aim to have the same results as I have already done with the combine tool in ArcGIS (see the attribute table).

Count = number of pixels
19911m = first raster (from 1951) with my typology code
2018 = second raster (from 2018) with my typology code

This methode help me to see transitions (eiher stable on unstable) between the two rasters, between two dates.

As exemple, on the attribute table we can see :

127 (typology code ; raster 1951) to 127 (typology code ; raster 1951) = 17207 pixels

210 (typology code ; raster 1951) to 310 (typology code ; raster 1951) = 108264 pixels

PS: all my rasters have the same typology applied, same resolution, same area.


texts: Are there any ways to composite Sentinel-2 bands with python?

texts:
Sentinel-2 satellite bands have different spacial resolutions, so how can I composite bands of satellite image in python (instead of gdal_merge script)?
 I tried to stack them by numpy, write each raster in each output band by GDAL, resize and reshape rasters and stack them. All of them are not working.
When I tried to stack bands in QGIS, it works perfectly, all bands are available.


texts: Difficulty using ordinary kriging in QGIS 3

texts:
I have a shapefile with 65 data points. Each data point has a value for NO2 pollution. I would like to use kriging to create a grid of NO2 pollution values  for the whole area. I have tried using Processing> Toolbox> SAGA> Ordinary Kriging. The process runs OK and two new layers are created called Quality Measure and Prediction, but both these layers are invisible or show no information.

What I assume is the value range of the Prediction file seems to indicate that something has gone wrong. It shows that the value of the black shade is 1.79769e+308 and the value of white is -1.79769e+308.

Can anyone help?


texts:
I am importing a KMZ we receive as output from an iOS app into ArcMap using the KML to layer tool.  The resulting table includes a "pop-up info" column that contians cells with marked up tags like:

Description	Stakes for monitoring 2 stakes

Time	N/A

Description	N/A

Is there a simple way to strip these html tags with the field calculator and VB Script / Python within ArcMap, or from the KMZ itself somehow prior to importing?

Thanks!

would a period at the end of a sentence interrupt the script?  That's the only thing I noticed between the stopping point cells...

Description	View northwest toward point between alignment and high bank. Standing on alignment on left (left) bank of stream channel.

Time	N/A

Photos


texts: Why does QGIS 2.4 prematurely terminates on startup at "Starting Python"?

texts:
After reinstalling Python (3.4) as well as OSGeo4W and the QGIS Standalone, my QGIS always crashed when "Starting Python" on my Windows 7 64 Bit environment. It simply did not proceed to the QGIS-GUI to be of any use. Furthermore, this happened on the way:

As I was beginning to write a plugin for QGIS, I faced a problem where the ui-compiler for qt-uis wouldn't work without telling it where the PYTHONHOME and PYTHONPATH is. As advised in a post here on gis.stackexchange, I set the respective system variables and added PYTHONHOME to my PATH-variable as well. The two steps ui-compilation worked, the test-plugin was finally ready but... QGIS refused to start as mentioned.

No errors are thrown/logged.

How do I get Quantum GIS 2.4 Chugiak back to work?


texts: Performance issue with SQLite + Spatialite + C#

texts:
I am having a performance issue with SQLite database (.db) which extends to SPATIALITE

I am trying to update Geometry of 1,00,000 records in database (.db) which taking around 50 minutes.  Too much slow.

my code is like below ::

        for (int q = 0; q < list.Count; q++)
            {
        ArrayList castarraylist = new ArrayList();
        castarraylist = (ArrayList)(list[q]);

        using (var cmd = new SQLiteCommand(con))

            using (var transaction = con.BeginTransaction())
            {
                cmd.Transaction = transaction;

                for (int y = 0; y < castarraylist.Count; y++)
                {
                        cmd.CommandText = Convert.ToString(castarraylist[y]);
                           cmd.ExecuteNonQuery();
                }
                transaction.Commit();
                GC.Collect();
            }
        }

Here each castarraylist contains 5000 records. which updating into database with transaction. so loop go through 20 times and complete the update all.
While I manually check the time it's increasing the time at each iteration for 5000 records. like

1st 5000 records processing time > 1:11 minute

2nd 5000 records processing time > 1:25 minute

3rd  5000 records processing time > 1:32 minute

4th 5000 records processing time > 1:40 minute

5th 5000 records processing time > 1:47 minute

6th 5000 records processing time > 1:52 minute

...

...

...

17th 5000 records processing time > 3:32 minute

18th 5000 records processing time > 3:44 minute

19th 5000 records processing time > 4:02 minute

20th 5000 records processing time> 4:56 minute

Why this happening I don't able to understand.
My sourcecode written in C# and my laptop configuration is i5 2.6 GHz, 4 GB RAM, 500 GB HD.

I made connection like below ::

SQLiteConnection con = new SQLiteConnection("Data Source=" + fullPath + ";Version=3;Count Changes=off;Journal Mode=off;Pooling=true;Cache Size=10000;Page Size=4096;Synchronous=off");

(*fullpath - is my database path)


texts:
I would like to read one variable (column) from a very large gpkg file (20GB) into R. It has only one layer. Is it any way to get one column (I assume that it is a character type but it could be a numeric one, I am not sure how this data is stored).
I tried this way:

df_pi <- st_read("All_DB2_Shapes.gpkg", query = "SELECT POLY_ID")

But I received this message:

  Error in CPL_read_ogr(dsn, layer, query, as.character(options), quiet,  :
   SQL execution failed, cannot open layer.
  In addition: Warning message:
  In CPL_read_ogr(dsn, layer, query, as.character(options), quiet,  :
   GDAL Error 1: In ExecuteSQL(): sqlite3_prepare_v2(SELECT POLY_ID):
   no such column: POLY_ID

If I just try to open it with st_read command then it cannot manage it because of the lack of memory.


texts: Growing error from west to east when overlaying raster and shape

texts:
I am using a shape of Austria from gadm.org which I want to overlay with a netCDF file in R.
As the netCDF has projected coordinates in EPSG:31287, I have reprojected the shape using FWtools and the .proj4-string from spatialreference.org.

Unfortunately the shape does not fit the raster data: The best fit is in the west, and the error grows in eastern direction (shape has a bigger extent in eastern direction than the raster). The north/south extent is correct.
Any ideas what could be wrong?
I already tried to play around with lon_0 values, but as far I am informed correctly, this causes only a general shift.


texts: Extrude between wrong on the edge, ArcScene

texts:
I'm trying to extrude between two Tins which are equally shaped. However if Extrude Between those tins using input feature created by minimal_bounding from one of them I get a strange error that two sides of my extruded block are inclined and thus the result is wrong (please have a look at the img attached).

I have absolutely no clue what could be the reason. I would be very glad for comments. Thank you.

Here are some images to illustrate the problem (thanks for the reputation so I can finally post them ;) ):

a http://imageshack.us/a/img833/4924/48164856.jpg
b http://imageshack.us/a/img801/706/94039199.jpg

I'm using ArcScene with ArcGIS 10.1.

This question can also be found on:
http://forums.arcgis.com/threads/77930-Extrude-between-wrong-on-the-edge-ArcScene?p=273925#post273925


texts: Creating buffers on selected point layer in ArcGIS Desktop?

texts:
I try to create buffers of 5m on the selected points only (in ArcGIS 10.3):

The problem is that the buffer tool create buffers for all the features,  why is that?


texts: Extracting features from LiDAR point cloud with ArcGIS: 'no spatial reference exists'?

texts:
I have LAS data obtained from 2014 NCFMP Lidar: Statewide North Carolina (Phase 2). I've created LAS dataset in ArcGIS and was trying to extract building footprints using LAS to Raster tool (I filtered classification codes first). However, when I ran the tool it displayed an error 'ERROR 999999: Error executing function. No spatial reference exists.' I've looked at metadata and horizontal coordinate system seems to be correct while vertical should be NAD83(2011) Ellipsoid (meters) and in properties of my LAS dataset it displays as Custom NAVD88 height. I reassigned to NAD83, but it didn't change anything and I'm not sure how to fix it otherwise.

I  also opened my data in ERDAS Imagine to double check, and, indeed, elevation information was shown as undefined.

I'm new to LiDAR.

I've attached screenshots in case it'll be helpful.

  Metadata from 2014 NCFMP Lidar:

  Process Description: Data were received by the NOAA Office for Coastal Management from the North Carolina Flood Mapping Program in the form described above. It was ingested into the Digital Coast Data Access Viewer system for distribution by:
    1) reprojecting to geographic coordinates with vertical meters
    2) transforming the vertical datum to NAD83 ellipsoid heights using the NGS GEOID12a model
    3) compressing the data using laszip (LAStools version 150828) to LAZ format.

Process Date: 201509
Direct Spatial Reference Method: Point

*Spatial Reference Information:
**Horizontal Coordinate System Definition:**

Geographic:
Latitude Resolution: 0.0000001
Longitude Resolution: 0.000001
Geographic Coordinate Units: Decimal degrees
Geodetic Model:
Horizontal Datum Name: NAD83 (2011)
Ellipsoid Name: Geodetic Reference System 80
Semi-major Axis: 6378137.000000
Denominator of Flattening Ratio: 298.257222

**Vertical Coordinate System Definition:**
Altitude System Definition:
Altitude Datum Name: NAD83(2011) Ellipsoid
Altitude Resolution: 0.001
Altitude Units: Meters
Altitude Encoding Method: Explicit elevation coordinate included with horizontal coordinates


texts: How to include layer names when union'ing in virtual layers?

texts:
I am creating a virtual layer in QGIS using the following statement:

SELECT * FROM (
SELECT * FROM a UNION ALL
SELECT * FROM b UNION ALL
SELECT * FROM c
) AS merged_tables

What is the best way to add a column to merged_tables containing the layer name of the source tables (in this case: a, b, c)? It would also be nice to be able to include the path to the source for a, b and c.


texts: Finding Percentiles from a List in Google Earth Engine

texts:
I am using Soil Moisture data from the FLDAS model, and have reduced the data down to a list containing a series of means. These means are average Soil Moisture for the Kenyan dry season (January through March 1982-present).

I wish to find the 20th, 30th, and 45th percentile value from this list using the ee.Reducer.percentile() function, but I keep getting the following error:

"ComputedObject (Error) List.reduce: Input must be a scalar number."

Looking for any help/ advice on finding percentile values from a list of numbers. Here is the code I have been using:

//importing collection and defining boundary of interest
var boundary =ee.Geometry.Point(37.559794, 0.771817).buffer(10000);
var SM = ee.ImageCollection('NASA/FLDAS/NOAH01/C/GL/M/V001').select("SoilMoi00_10cm_tavg");

//years of interest
var years = ee.List.sequence(1982,2019);

//reducing to Short Dry season
var temp = ee.ImageCollection.fromImages(
  years.map(function(y){
      return SM
        .filter(ee.Filter.calendarRange(y, y, 'year'))
        .filter(ee.Filter.calendarRange(1, 3, 'month'))
        .mean()
        .set('Season', "Short Dry").set('year', y)
        ;
  }).flatten());
  print("Check Star/End Dates", temp);

//Finding mean pixel value of each yearly image (reduces to single number)
  var temp = temp.toBands();
  var temp = temp.reduceRegions({
      collection: boundary,  // the region over which values are sumamrized
      reducer: ee.Reducer.mean(),  // the summary statistic
      scale:1000 }).first();

//Putting these numbers in a list
  var propNames = temp.propertyNames();
var list = propNames.map(function(image){
  return temp.get(image);
});

//Calculating 20th,30th,and 45th percentiles
var percentiles = list.reduce(ee.Reducer.percentile([20,30,45]));

print("Percentiles",percentiles);


texts: Using os.remove in ArcGIS Python script tool?

texts:
I have a Python script that I want to turn into a scripting tool in ArcGIS.  Part of the script creates several SDE database connection files, creates a versions using those connection files, and then deletes the connection files.

To delete the connection files, I'm using os.remove(connections).  This works as expected in PyScripter and the Python window in ArcMap.  However, when I create a scripting tool and try to run it, it returns the message, "System cannot find the path specified" in reference to this line.

I can modify the code to say

if os.path.exists(connections):
    os.remove(connections)

so that it will move on to the next step, but then  I am left with a folder full of connection files when the tool is finished running.

Any suggestions?

I tried arcpy's Delete, but got Error 000601 that it may be locked by another application.  I found this bug which says that this is a bug with arcpy.Delete_management in ArcMap 10.2.1 (which I am using), and suggests using os.remove() instead.


texts: Return field as array from any map layer in QGIS expressions

texts:
In QGIS 3.8 It's obvious in the expression builder how to reference a layer or how to get a feature from a layer but not how to return a whole field from a layer.

What expression will return a field of a map layer as an array?


texts: Merging .dbf Tables in Python w/o merging all fields of same name

texts:
I'd like to merge several tables together using Python. They all have the same field names: "TIME", "X", "Y", and "VALUE".

I would like the "VALUE" fields to remain separate when merged, so that the new table has the fields: "TIME" "X" "Y" "VALUE1" "VALUE 2" "VALUE3"...and so on.

I wrote a script based off the ERSI script for merging, and tried to use Field Mapping in a way that would let me rename "VALUE" for each table before merging.

# Name: Merge.py
import arcpy
from arcpy import env

# Set environment settings
env.workspace = "C:/Users/user1/Documents/Processing"

# Get tables to be merged
listTable = arcpy.ListTables("*", "ALL")

# Create FieldMappings object to manage merge output fields
fieldMappings = arcpy.FieldMappings()
i = 1

# All tables have a column "VALUE" only I don't want them to merge
for tables in listTable:
    fieldMappings.addTable(tables)
    fldMap_SOSvalue = arcpy.FieldMap()
    fldMap_SOSvalue.addInputField(tables,"VALUE")
    SOSvalue = fldMap_SOSvalue.outputField
    SOSvalue.name = "SOSvalue" + str(i)
    fldMap_SOSvalue.outputField = SOSvalue
    fieldMappings.addFieldMap(fldMap_SOSvalue)
    i += 1

# Use Merge tool to move features into single dataset
mergetable = "C:/Users/breezevi/Documents/ConModelTest/Tables/2013.dbf"
arcpy.Merge_management(listTable, mergetable, fieldMappings)
print "Merge is done"

This script runs, however, it does not seem to merge correctly. The resulting table has the fields "TIME" "X" "Y" "VALUE1" "VALUE2" "VALUE3" ... but all the numbers in each "VALUE#" column seem to be the same (e.g. it reads 1 | 1 | 1 | 1 across the whole row, but it should read 1 | 3 | 5 | 2 ). Almost like the script is creating all the new fields, but only filling them with values from one table. Only, I'm not sure where I went wrong in the script.


texts: Labelling layer not from attribute fields QGIS API v2.14?

texts: How to label vector layer with my own values that aren't stored in the attribute fields?


texts: Distinguish between point and any other spatial type

texts:
When running an ST_INTERSECTION query on two linestrings in PostGIS, I expect it to return a point (e.g. 'POINT(1 2)'). This works 99,99 % of the time.  However, sometimes (for whatever reason), my query returns e.g. a MultiPoint, LineString, ... or any other spatial type.  In that case, I simply want to skip my queryResult and move to the next row in my table.

So I am looking for something like :

  if queryResult.SpatialType() = POINT:
         'do something'
  elif queryResult.SpatialType() != POINT:
         pass

But I cannot get anything like this to work.

What would be the correct statement to make a proper distinction between a point as the queryResult and any other spatial type ?


texts: Understanding the output of "Table Compare" in ArcMap

texts:
I'm running ArcGIS Desktop 10.5.1 on Windows 10. I tried using the "Table Compare" tool to compare the attribute table of an old version of a shapefile with a the attribute table of a newer version of the same shapefile. The newer version contains some added entries (old 1898 rows / new 2088 rows), this is the only difference I can tell when I manually compare the tables.
When checking the output file I got something like this:
Has_error, Identifier, Message, Base_value, Test_value, ObjectID
"true", "Table", "Table row counts are different", "1898", "2088", -1
"true", "Table", "Table row counts are different", "1898", "2088", -1
"true", "FeatureClass", "ObjectID 0 is different for Field Shape", "Geometry", "Geometry", 0
"true", "FeatureClass", "ObjectID 1 is different for Field Shape", "Geometry", "Geometry", 1
"true", "FeatureClass", "ObjectID 2 is different for Field Shape", "Geometry", "Geometry", 2
"true", "FeatureClass", "ObjectID 3 is different for Field Shape", "Geometry", "Geometry", 3
.
.
.
"true", "FeatureClass", "ObjectID 1897 is different for Field Shape", "Geometry", "Geometry", 1897
"false", "Table", "Tables have same number of fields", " ", " ", -1
"false", "GeometryDef", "GeometryDefs are the same", " ", " ", -1
"false", "Field", "Field properties are the same", " ", " ", -1

The first two (three) lines make sense to me, even if I cannot understand why this is listed twice. The following 1898 lines regarding the field "Shape", however, make no sense to me.
EDIT: When I check the attribute tables of the two shapes in ArcMap, all the entries in the field "Shape" are "Polyline", not "Geometry".


texts: Sending a https request from a QGIS plugin using urllib2

texts:
I'm developing a QGIS plugin that needs to send https request to an API.
I have to use urllib2 (see below if you want to know why I don't just use the requests package)

The first request I need to make is a POST one, with id in the header, so the API would grant me an access token (according to the OAuth2 protocol)
I have this portion of script, that works when executed alone, but doesn't as soon as I integrate it in my QGIS plugin script, it gives me a weird URLError (code 10054 : the distant host had to close a connexion) but I suspect QGIS doesn't event send the request at all, for it isn't an answer that the API gives.

import urllib
import urllib2

myIDs = "Basic hjNjkJCOHdrUncN4RVV"
url = 'https://id.mygreatapi.com/oauth/token'
values = {"grant_type":"client_credentials"}
data = urllib.urlencode(values)
header = {"Authorization": myIDs}
req = urllib2.Request(url,data, header)
try:
    handle = urllib2.urlopen(req)
    content = handle.read()
    print(content)
except urllib2.HTTPError as e:
    code = e.code
    texte = e.read()
    print(code)
    print(texte)

Note that when I copy a valid token into my QGIS plugin script and then make GET requests to the API using http, it works. So I guess the problem comes from the https.

Any ideas why QGIS wouldn't send this request ?

Why don't I use requests library ?
It is necessary that the plugin could be installed by anyone (especially somebody without admin rights on his computer, and no programming knowledge whatsoever) so I can't expect people to use their OSGeo and use pip install.
The solutions proposed in this page being totally out of my competences, I considered that I have to stick with the packages that are built-in the python version QGIS uses.


texts: Extend the contours of the polygon into a polygon

texts:
I would like to be able to extend the vertices of one polygon to the edge of another polygon, how does the snapping between poligons using python.
The photo shows a representation of my problem.

I would like the orange polygon snap the red contour (the other polygon), as in the photo below.

I am using fiona as shown in this example (How to do snapping in Geopandas?),  but I couldn't get the result.
from shapely.ops import snap
import geopandas as gpd

border_shp = gpd.read_file(...)
multipolygon = gpd.read_file(...)

for index, row in multipolygon.iterrows():
    closest_geom = list(border_shp ['geometry'])[1]
    snapped_geom = gpd.GeoDataFrame(snap(row['geometry'], closest_geom, 8))

I would be grateful for the help


texts: Is it possible to see the selected geometry in the map composer? And change the appearance of them?

texts:
I have QGIS 3.4.4 Madeira. I would like the selected geometry to be seen in the map composer, but I can't find the way to do it.
Also, I would like to change the appearance of the selected geometry, with a outline border and no fill inside. I know how to change the fill color of selected items and how to set transparency, but this is not what I need.
In ArcMap is quite easy to do, but I can't find the way in QGIS.

EDIT:
This image is what I want to do in QGIS, but have been done in ArcMap. The selected plot have a symbology in red lines with transparent bottom (ArcMap allows to set a different symbology to selected geometry in each layer). The green color indicates the urban classification of this area. The rest of the plots don't appear because they have no symbology. I want to have this on the map composer with the selected plot or plots, print the map to PDF, select another plot, print the map to PDF, an so on. The plot to select don't attend to any rules, just to people asking information about them.


texts: Retrieving values from grid squares using QGIS?

texts:
I have a map of Europe, with grid squares that each have vales for temperature. I need to have a list of all the temperature values for each country.
Instead of hovering my mouse over each relevant grid square and manually writing down the values, is there a way I can select and retrieve a list of all in the selection?

It's also ECMWF data if that helps.


texts: Point for Each Parcel on the Road in front with angle information

texts:
I have this interesting problem and hope someone can provide some insight/solution. I have access to full blown ESRI software.

I have property parcels layer with close to 20,000 records. The attributes of this layer contains regular parcel info, plus addresses. There are three separate columns for addressing; Street Number, Street Name, and Full Address.

I want to place a point right at the centre front of each parcel, preferably on the Road line feature that I have. In case where the the parcel is on the corner and has access to more than one road, I want to place the point on the road that matches with the address. Finally for each point I need an xy coordinate and an angle for facing the front of the parcel.

This data will then be used by a developer to incorporate Google Street View with a property information application.

Does anyone have some ideas as to what would be the best way to accomplish this?


texts: Place all labels on a map non overlapped

texts: As asked for ArcMAP by @dassouki, is there a way to place all labels on a map non overlapping using QGIS?


texts: Determining distance between two point locations specified in degrees, minutes and decimal seconds?

texts:
Without using GIS software, how can I determine how far it is from 20°01'12.07" N -75°47'53.57" W to 20°01'15.28" N -75°47'46.02" W, which are two point locations (historical battles) specified in degrees, minutes and decimal seconds?

Point A: 20:01:12.07N; 75:47:53.57W
Point B: 20:01:15.28N; 75:47:46.02W


texts: Heat mapping forestry height trends using ArcGIS?

texts:
I am using ArcGIS and have access to Spatial and 3D Analyst extensions.

I've attached a link to a picture of a 1m raster “canopy height model” which was derived from LiDAR and displays tree heights (in 5m intervals but can  obviously be adjusted). I do 3d models w/ a program that uses polygon inputs to derive height classes.

How can I convert this raster to a polygon and capture broad polygon classes?

(I can’t use raster to polygon as I work with 100’s of hectares at a time and an individual polygon for each raster cell would be far too much detail and slow the program down.)

I’m thinking of a heat map (?) type exercise which captures general trends as polygons. I’ve circled several polygons in blue that would be the ideal type of polygons I would like to capture…but ultimately I do need 100% polygon coverage.

I have access to SA and 3D analyst if that helps! In terms of accuracy, I'm not sure. As accurate as I can be without having 1m polygons or anything too small.


texts: Splitting map, island tied in same polygon as mainland using QGIS

texts:
I have a map originally as a shapefile. The map includes a country (mainland) and an island far away from the mainland. I want to cut the map so that only the mainland is shown and save it again as a new shapefile. The problem is that the island is tied in the same polygon as one of the mainland regions.

How would you advise me to split the map?

I am using QGIS 3.10.5.


texts: CARTO update existing table from external MySQL database

texts:
I would like my Carto data to be updated, using external data (MySQL database, I run and manage locally using PHP scripts), so that the Carto table gets updated with the most recent data, whenever the map gets loaded.

Is it possible?


texts:
When I try to connect to a small localhost PostGIS database using QGIS, I get the error message:

FATAL: no pg_hba.conf for host "fe80::4528:3236:418a9e7b%10" user "postgres", database "nyc", SSL off

The file pg_hba.conf contains

# TYPE  DATABASE        USER            ADDRESS                 METHOD

# IPv4 local connections:
host    all             all             127.0.0.1/32            md5
# IPv6 local connections:
host    all             all             ::1/128                 md5
# Allow replication connections from localhost, by a user with the
# replication privilege.
#host    replication     postgres        127.0.0.1/32            md5
#host    replication     postgres        ::1/128                 md5

I can connect just fine using pgAdmin. I noticed that a few others have asked similar questions here, but I can´t find a solution using localhost. Any suggestions on what to try? QGIS 1.8.0 and PostgreSQL 9.2/PostGIS 2 on Win Vista 64 bit.


texts: Loading points and polylines into PostGIS database and creating connectivity between them using topology?

texts:
I am new in PostGIS, and I am trying to create connectivity(using topology) between 3 sets of points and polyline to model an electrical distribution network.

Does anyone know how to do this?


texts: Proj4js Problem in converting coordinates from google mercator

texts:
I am a really beginner into web-cartography, so sorry if this subject is rudimentary or basic but I truly search for a long time and I am still stick in it.

So the problem is:

I don't understand why this code don't convert correctly my coordinates (exprimed in google mercator epsg:900913) to WGS84 web mercator epsg:4326:

   Proj4js.defs["EPSG:4326"] = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";
Proj4js.defs["EPSG:900913"]= "+title=GoogleMercator +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs";

var dest = new Proj4js.Proj("EPSG:4326");
var source = new Proj4js.Proj("EPSG:900913");
var point = new Proj4js.Point( 85286.417057415703, 4347668.9422677439 );
Proj4js.transform(source, dest, point);
var transformation = Proj4js.transform(source, dest, point);
alert(point.y + "," + point.x);

I should have in return x (long) = 0.7661409 and y(lat) =36.3411771 (calculated by http://twcc.free.fr/) but with this code I have:
x=0.0003264583485814848 y=0.000006882360979601046


texts: I'm trying to simplify a polygon using the interior angle of a vertex as the "tolerance" instead of the distance between vertexes. E.g. the polygon in the image below would lose one vertex, because of its large internal angle. As an example, I might want to remove all vertexes with an internal angle greater than 170°.


texts: QGIS Field Calculator Update Existing Field

texts:
I am using QGIS 2.18 and am trying to replace the values in an existing column. The errors I get are that Replace doesn't have the right number of arguments and the no root node error.

This is likely a python thing and I'm not that experienced in python. The example indicates three parts: input string, replace string, the output string.  I want to replace the whole thing.


texts: Aggregating point data regarding lat/lon and year

texts:
I have two datasets:

The first one contains georeferenced point data. The structure looks like that:

The second one contains georeferenced point data + radius. The structure looks like that:

My goal is to get the sum of the "newusdco" for the points in the corresponding area of the second dataset.

My first problem is, that I don't know how to aggregate the data to the lat/long+radius data.

My second problem is, however, that the data varies over time too. Meaning that I only want to get the sum of "newusdco" in the conflict site in the same years.

I have been reading up in the forum a lot but couldn't find a solution, since I am not very familiar with R.


texts: Atmospheric Calibration for Landsat 7 ETM+ Gap Filling

texts:
This is for thesis work involving mapping of above-ground mangrove forest biomass using Landsat 7 ETM+ imagery.

I am assessing correlations between biomass with vegetative indices, e.g. NDVI, and then mapping biomass. I am using a 04/02/2013 Landsat 7 ETM+ image with Scan Line Corrector off, so I need to gap fill. I am using a 08/03/2013 image for gap fill.

Both images have been atmospherically calibrated before gap fill using DOS_improve. However, the reflectance of the two images still contrast to each other. I do not know why the reflectance is contrasting, or how to better match the reflectance of the 08/03/2013 image to the 04/02/2013 image.


texts: Cloud and Haze Filter Not Working for Landsat 5 Data in GEE

texts:
I'm trying to make an annual composite of Landsat 5 data for a region of interest in Google Earth Engine but what I end up with is hazy and cloudy. I've tried filtering for haze and clouds but still get both in my final version. I'm so confused as to why this isn't working and would really appreciate any guidance that anyone can offer me. The code I have written is below

 roi =
ee.Geometry.Polygon(
    [[[-70.94315884724863, -13.228382487713638],
      [-68.53715298787363, -13.228382487713638],
      [-68.53715298787363, -12.022330477798716],
      [-70.94315884724863, -12.022330477798716],
      [-70.94315884724863, -13.228382487713638]]]);
var clouds = function(image) {
 var qa = image.select('pixel_qa');
 var cloud = qa.bitwiseAnd(1 << 5)
              .and(qa.bitwiseAnd(1 << 7))
              .or(qa.bitwiseAnd(1 << 3));
 var mask2 = image.mask().reduce(ee.Reducer.min());
 return image.updateMask(cloud.not()).updateMask(mask2);
};

function lsHazeMask(image) {
var mask = image.select('sr_atmos_opacity').lte(300);
return image.updateMask(mask);
}

var jsonCoordString = roi.toGeoJSON();

for (var start_year = 1986; start_year <= 1990; start_year++)
{
var startDate = ee.Date.fromYMD(start_year,01,01);
var sr = ee.ImageCollection(ls5SR).filter(ee.Filter.bounds(roi)).filter(ee.Filter.calendarRange(6,9,'month')).filterMetadata('CLOUD_COVER', 'less_than', 30).filter(ee.Filter.date(startDate,startDate.advance(1,'year')));
sr = sr.map(clouds);
sr = sr.map(lsHazeMask);
sr = sr.reduce(ee.Reducer.median());

Map.addLayer(sr.select([2,1,0]).clip(roi),{min: [0,0,0], max: [2000,2000,2000]},'puerto'+start_year);
}


texts: How can I edit spatial data automatically on postgis and check for topology errors on said data whitout using shapefiles?

texts:
So, I have spatial data stored in postgis (points, lines and polygons) and would like to edit and update those in an Open Source software (quantum GIS, grass, etc.) to check for topologic errors, among others stuff, (arc-node-polygon topology) in an automatic way. So when I edit or update this data it would automatically be saved in the geodatabase whithout using the
import from DB -> edit -> create shapefile -> import shapefile to DB routine.
I know the ArcSDE can do all of this, but I need it done with Open Source software alone. Is there a way to do it?

Ps: It would be best if those editions and updates would respect my database relationships between tables, but if not, the problem can be solved using triggers on postgis.

Thanks in advance!


texts: Finding spatial data on Enumeration Area boundaries from South Africa 2001 Census?

texts: I am trying to locate spatial data for the boundaries of South African Enumeration Areas (EAs) from the 2001 Census. The links to the data from the 2001 census on the Statistics South Africa website say "object not found".


texts: QGIS Python Console AddMaplayer but with low resolution (very dim)

texts:
Update： Figured. I used the wrong style file.
I'm new to QGIS Python develop. I was trying to add raster layer (.asc file) into QGIS. Below are my srcrips.
from qgis.PyQt.QtCore import *
import gdal
import os

filePath = "path/2017011003.asc"
style_layer = "path/StreamStyle.qml"

fileInfo = QFileInfo(filePath)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(filePath, baseName)
if not rlayer.isValid():
  print "Layer failed to load!"

crs = rlayer.crs()
crs.createFromId(4610,QgsCoordinateReferenceSystem.EpsgCrsId)
rlayer.setCrs(crs)

rlayer.loadNamedStyle(style_layer)
QgsMapLayerRegistry.instance().addMapLayer(rlayer)

I did add the raster layer into qgis, but the problem is that the layer seems very dim... I also import the same raster image by just using the qgis "Add Raster Layer" function which shows quite clear
Here is the comparison:
added by scripts (dim)

added by using QGIS "Add Raster Layer" function


texts:
I am totally new to doing this using QGIS. Can anyone help me to produce the following, which I believe is commonly known as a ZTV analysis.

I am running Wein and have the Viewshed Analysis plugin loaded. I have downloaded the asc data for a single UK grid square, extracted it, converted it and then saved it as a Geotiff.

This is what I have been asked for:

We used the free 50 mm Digital Terrain Model (DTM) from the OS, gave the proposed development a height e.g. approximately 8 metres for a 2 storey house within the site boundary and viewing height of approximately 2 metres. We would sometimes ‘screen’ the ZTV by introducing shape files for areas of woodland and settlement

I really do not have a clue where to start.


texts:
I am trying to load a vector file to QGIS Python console.

v_layer = iface.addVectorLayer("G://Consultancy\Data\Transfer\Bodhi (Louis) Gorringe\QGIS\Shapefiles\Biological Shapefiles\Seabirds\At sea densities of gannet in the breeding season\710_br.shp")

Below identifies the error code I always get when I enter this script.

From looking at other examples online it seems that they add a bit at the end of the script within single quote marks ('). As seen in the image below in the red box.

I am aware this is specified for where they have saved their documents, but I am unaware of what the (') columns mean. I feel this is the reason why my vector layer wont load. Can anyone confirm, or know a reason why the vector layer wont loads?


texts: Create Choropleth Style Map in CartoDB programmatically

texts:
I am trying to create a choropleth map for a school project I am working on. Based off user input, I will create a csv file that I will then upload to CartoDB using their Import API. By setting the create_vis=true option, CartoDB will automatically make a map using that table. By default, however, this map is of the Simple type. Is there any way to make this a choropleth type map programatically?

I understand that you can do this in the wizards editor in the UI, but I am trying to do all of this at runtime (so that it is dynamic).


texts:
This script should make a field for all the feature classes in a geodatabase, populate the new field with the name of the feature class it's being created in, and then merge all the feature classes together. However, when I run the code I get the success message but only the first feature class features are in the output of the merged feature class.

import arcpy, os

arcpy.env.workspace = r"J:\Non_RS_MXDs\Crown_Scrap.gdb"

for fc in arcpy.ListFeatureClasses():
    arcpy.AddField_management(fc, "DyRt", "TEXT", field_length = 30)
    with arcpy.da.UpdateCursor(fc, "DyRt") as cursor:
        for row in cursor:
            row[0] = fc
            cursor.updateRow(row)

fClasses = []
i = 0

try:
    fcs = arcpy.ListFeatureClasses()
    for fc in fcs:
        print fc
        fClasses.append(fc)
        i = i+1

except Exception as err:
    arcpy.AddError(err)
    print err

try:
    arcpy.Merge_management([fClasses], "Rts_8_15")
    print "merge: SUCCESS!"

    #count number of new features
    result1 = int(arcpy.GetCount_management("Rts_8_15").getOutput(0))
    print result1

except Exception as err:
    arcpy.AddError(err)
    print "merge: FAIL :("
    print err


texts: Using dynamic property value in OnlineResource href against MongoDB store?

texts:
I am passing this SLD to geoserver:

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd">
<NamedLayer>
    <Name>cite:test</Name>
    <UserStyle>
      <Title>Test</Title>
        <FeatureTypeStyle>
          <Rule>
            <PointSymbolizer>
              <Graphic>
                <ExternalGraphic>
                  <OnlineResource xlink:type="simple" xlink:href="http://localhost/icons/river"/>
                    <Format>image/png</Format>
                </ExternalGraphic>
              <Size>35</Size>
            </Graphic>
          </PointSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

That works great, but will get the my 'river' icon for ever feature.  Say there is a property called icon.type in for the data source that lives behind this layer.  Is there a way to use the value of that property for each features OnlineResource?  In that way I would only have to define one PointSymbolizer.

IE, something like this:

<ExternalGraphic>
  <OnlineResource xlink:type="simple" xlink:href="http://localhost/icons/{icon.type}"/>
  <Format>image/png</Format>
</ExternalGraphic>

EDIT:

Turns out you can indeed use a dynamic property by using CQL in the url:

<ExternalGraphic>
  <OnlineResource xlink:type="simple" xlink:href="http://localhost/icons/${icon.type}"/>
  <Format>image/png</Format>
</ExternalGraphic>

However I am using a mongodb store and the extension does not seem to query the backend mongo store correctly for that property.


texts: How to move polygon pie charts symbols in ArcGIS?

texts:
I used pie charts in my map, but they don't appear in the desired position. Instead, they all appear in the exact middle of polygons.

The charts were created via symbology -> charts -> pie option

I need to put area names to the middle of the polygons and the pies above them, or adjust the positions of both of the text and pies. But I can't seem to figure out how to move the charts, not even with Editor.


texts: Move/Shift a line toward x and y axises in SQL-Server

texts: Is there a function in SQL Server to move/shift a geometry toward X or Y axises?


texts: Function-based spatial index: Use sdo_geometry() constructor instead of custom function?

texts:
I have an Oracle 18c table that has XY columns, but not a SHAPE column. I've created a query on the table that dynamically generates an sdo_geometry SHAPE column.
sdo_geometry(2001, 26917, sdo_point_type(longitudex,latitudey, null), null, null) as shape

Function-based spatial index:
In an effort to improve the performance of the query, I've created a function-based spatial index.
The Oracle docs provide instructions for creating a custom function for the purpose of the function-based index. So that's what I've done:
--Create a deterministic function:
create or replace function sdo_geom_from_xy(longitudex in number, latitudey in number)
return sdo_geometry deterministic is
begin
    return sdo_geometry(2001, 26917, sdo_point_type(nvl(longitudex,0), nvl(latitudey,0), null),null, null);
end;

--Update USER_SDO_GEOM_METADATA:
delete from user_sdo_geom_metadata where table_name = 'WO_MV2_SBX';
insert into user_sdo_geom_metadata values(
    'WO_MV2_SBX',
    'maximo.sdo_geom_from_xy(longitudex,latitudey)',
    sdo_dim_array(
        sdo_dim_element('X', 663795.93, 676674.94, 0.001),
        sdo_dim_element('Y', 4857659, 4879068.83, 0.001)), 26917);

--Create a function-based spatial index:
drop index wo_mv2_sbx_idx;
create index wo_mv2_sbx_idx on
    wo_mv2_sbx(maximo.sdo_geom_from_xy(longitudex,latitudey))
    indextype is mdsys.spatial_index;

--Create a view that uses the the function-based spatial index:
create or replace view wo_mv2_view_sbx as (
select
    objectid,
    maximo.sdo_geom_from_xy(longitudex,latitudey) as shape
from
    maximo.wo_mv2_sbx
where
    plussisgis = 1
);
--Note: In order for the view get the benefit the function based index, the view *MUST* use the exact same syntax that was used in the function-based index.

Question:
When I use the custom function, everything works as expected. However, it seems strange to me that I would need to create a custom function, when the out-of-box sdo_geometry() constructor would do what I want.
Is there a way to use the sdo_geometry() constructor in a function-based index instead of custom function?
Edit:
There is a related post here:
Can I create an Oracle Spatial View from a non-spatial table? (with an excellent answer from @AlberGodfrind.


texts: How does ArcGIS compute the distance between two points with a non-equidistant projection?

texts:
This is a follow-up question to my previous one, Can you suggest some well-written introductory texts about coordinate system projections?

Let's assume I'm working with the CH1903 map projection, which for all I know is conformal, but not equidistant. Meaning, angles (shape) have been preserved, but not areas, distances, or scale. (At least these have not been preserved exactly). So far so good.

I'm wondering what kind of calculation ArcGIS performs when I now want to calculate the distance between two points. In ArcObjects, I could use the IProximityOperator interface as follows:

IPoint a = ...,
       b = ...;

double distance = ((IProximityOperator)a).ReturnDistance(b);

Question: When I'm working with a reference system that does not accurately preserve distances, what would ArcGIS do when I query it for the distance between two points (as shown above)?

Does it simply do some Pythagorean maths (a2 + b2 = c2) to get the distance, meaning the returned distance will only be as accurate as the projection allows?
Or will it do something more complicated, like some form of re-projection, to get a more accurate distance?

(The same question, but more generally: Once that geometries have been projected, does ArcGIS perform all calculations simply in Euclidean space, or does the used map projection still influence calculations of distances, angles, areas, etc.?)


texts:
This is the first time that I have worked with geospatial data and QGIS, so please be kind.

I was given a shape file that contains house districts for the state of FL.  This link provides instructions (for ArcMap users) to import the house districts into Tableau, which is my goal.

I was able to generate a file that contains the district, and coordinates, but I need the order column, and am not sure how to do this using QGIS (assumedly using an add-in).

"shapeid","x","y"
"0","-87.634896","30.866261"
"0","-87.634871","30.866344"
"0","-87.634851","30.866536"

The format that I am looking for is ...

"shapeid","x","y","order"
"0","-87.634896","30.866261","1"
"0","-87.634871","30.866344","2"
"0","-87.634851","30.866536","3"


texts: OSError: OSError: [WinError 126] The specified module could not be found - osm2geojson for Python

texts:
I installed osm2geojson with Anaconda following these instructions:https://pypi.org/project/osm2geojson/
Now, every time when I try to import osm2geojson into my Jupyter notebook I receive the following error message:
OSError: [WinError 126] The specified module could not be found
All my attempts to find a working solution failed.


texts: problem with viewing ASCII data files in ArcGIS

texts:
I have some problem viewing an ASCII file in a structured way in arcGIS. I have tried running this code but it seems problematic at least with me. I am unable to run this. I think the code needs modification. Can you please tell me the requirements to run this code... and some other help.
the link contains the file as well as the IDL source file
https://drive.google.com/?authuser=0#folders/0B8iAKVklcaMUVU5qY2NNVFBRbDQ


texts: I'm using QGIS 3.23-Bonn on a Win 7 machine. I'm trying to add the "inch" symbol to a label, but can't figure it out. I've tried using (Diameter) || '\' " but it doesn't work. However for length in feet this works (Shape_Length,0)||'\" || '\n' || Diameter. This displays two labels, stacked. One is in feet and one is in inches. All I need for another one is "inches".


texts: Does the ArcGIS Flex API simplify polygons by default?

texts:
I have a feature access map service that I use on an ArcGIS Flex web application.
Some polygons published by the feature service, seem to be "simplified" by the Flex API before they are drawn on the map.
The original shape of the polygons is ellipsoidal, like potatoes. But they are drawn like hexagons or octagons with summits and sharp angles.

How can I prevent this behaviour?


texts: Accuracy assessment in R - Calculation of User Accuracy

texts:
I want to assess the accuracy of the land cover classes classified by supervised classification of Landsat images. As reference data I use the aerial photography.

I sampled the classified Landsat data and at the same validation point, I identified the land cover class from aerial photography. I.e., each of my verification points have two attributes: one from Landsat (landsat), second from the aerial photography (reference).

I want to calculate an error matrix (contingency table), and accuracy assessment (overall accuracy, user's producer's accuracy) from the set of verification points.

I found a package greenbrown and asbio to assess the accuracy of classification. My final results of user's and producer's accuracy are howver switched between two packages.

Please, which package correctly calculate the values for user and producer accuracy?

Reproductible example

library(asbio)

# create dummy data
landsat <- c(1, 1, 1, 2, 2, 2, 3, 4, 5, 5, 5, 1, 1, 1, 2, 2, 2, 3, 4, 5, 5, 3, 3, 2, 2)
reference <- c(1, 2, 1, 2, 2, 2, 3, 4, 2, 2, 5, 1, 2, 2, 2, 1, 2, 3, 4, 5, 5, 3, 3, 2, 2)

# calculate Kappa statistics
asbio::Kappa(landsat,reference)  # Kappa(class1, reference)

# check out Kappa results
$ttl_agreement
[1] 76

$user_accuracy
1     2     3     4     5
75.0  58.3 100.0 100.0 100.0

$producer_accuracy
1     2     3     4     5
50.0  87.5 100.0 100.0  60.0

$khat
[1] 68.1

$table
       reference
class1  1 2 3 4 5
      1 3 3 0 0 0
      2 1 7 0 0 0
      3 0 0 4 0 0
      4 0 0 0 2 0
      5 0 2 0 0 3

# ----------------------------------------------------------------------
# make the same calculation with the greenbrown package
# ----------------------------------------------------------------------

library(greenbrown)
library(strucchange)
library(raster)
library(Kendall)
library(plyr)
library(bfast)
library(zoo)

# calculate the contingency table
tab <- table(landsat, reference)   data

# let's see the tab
tab

           reference
landsat   1 2 3 4 5
        1 3 3 0 0 0
        2 1 7 0 0 0
        3 0 0 4 0 0
        4 0 0 0 2 0
        5 0 2 0 0 3

# calculate the accuracy assessement
greenbrown::AccuracyAssessment(tab)

                    1        2   3   4   5 Sum UserAccuracy
1                 3  3.00000   0   0   0   6         50.0
2                 1  7.00000   0   0   0   8         87.5
3                 0  0.00000   4   0   0   4        100.0
4                 0  0.00000   0   2   0   2        100.0
5                 0  2.00000   0   0   3   5         60.0
Sum               4 12.00000   4   2   3  25           NA
ProducerAccuracy 75 58.33333 100 100 100  NA         76.0

The user and producer's accuracy are switched between two packages !
Please, which calculation and estimation of the user's and producer's accuracy is correct?


texts: Can you disable panning in ArcGIS story map?

texts: I have a story map set up for a specific small area and I want to keep people from being able to accidentally swipe and end up across the world. Is that possible?


texts:
I have a similar issue that reported here: R lidR does not read the CRS of my las file. lidR cannot read the CRS while lasinfo from lastools can.
readLAS("USGS_LPC_MS_Tupelo_UTM15_2015_15RYQ605365_LAS_2016.laz")
#> proj_create: crs not found
#> proj_create: crs not found
#> class        : LAS (v1.2 format 1)
#> memory       : 253.4 Mb
#> extent       : 761421.7, 762000, 3436500, 3438000 (xmin, xmax, ymin, ymax)
#> coord. ref.  : NA
#> area         : 0.84 kunits²
#> points       : 3.32 million points
#> density      : 3.97 points/units²

But the LAS file is not 1.4 (1.2) and the use of rlas 1.4 and lidR 3.2 from github don't help. The laz file can be found here


texts: How do you read/write a FeatureCollection from/to a json file?

texts:
I am trying to read in a json FeatureCollection (geojson 1.3.1) Python (2.7), add some properties, and write it back out to a json file (or string).

I can read the file, and modify the properties, but I end up with 2 nested FeatureCollection objects:

{
  "type": "FeatureCollection",
  "features": {
    "type": "FeatureCollection",
    "features": [
      {
        "geometry": {
          "type": "Polygon",
          "coordinates": [
            [

            ]
          ]
        },
        "type": "Feature",
        "properties": {
          "gval": 1.0,
          "hectares": 0.218681933167,
          "xid": 0,
          "id": 1.0,
          "area": 2186.81933167
        }
      }
    ]
  }
}

And I think I've worked out that I creating my FeatureCollection incorrectly:

path=r'd:\temp\test.json'

with open(path) as data_file:
    data = json.load(data_file)

print data

feature_collection = FeatureCollection(data)

print feature_collection

The feature_collection ends  with nested features, which is probably why the resulting geojson.dumps has the nested FeatureCollection
{"features": {"features": [{"geometry":

I can think of a number of ways to hack it, but what is the correct way to create that FeatureCollection? Modify feature attribtues, and then write it back out to a file?

Content of the test.json file (minus coordinate values):

{
  "type": "FeatureCollection",
  "features": [
    {
      "geometry": {
        "type": "Polygon",
        "coordinates": [

        ]
      },
      "type": "Feature",
      "properties": {
        "area": 2186.81933167,
        "id": 1.0,
        "xid": 0,
        "gval": 1.0
      }
    }
  ]
}

Version of geojson:

Metadata-Version: 2.0
Name: geojson
Version: 1.3.1
Summary: Python bindings and utilities for GeoJSON
Home-page: https://github.com/frewsxcv/python-geojson
Author: Corey Farwell
Author-email: coreyf@rwell.org
License: BSD
Location: c:\python27\arcgis10.2\lib\site-packages
Requires: setuptools


texts: Connect multiple 2 points with a line in ArcGIS?

texts:
I have a table in which each row has two columns representing the objectID of two different point features in the Map. There are thousands of such rows. How do I connect the each of these two points with a line?

I don't want the points of two different rows connected to each other.


texts: Updating all bands of PostGIS raster in one function (like st_setvalues) with one geomval[] set?

texts:
I have ~700 raster tables in postgres. Some have an expanse the size of California. All the tables are tiled and share tile IDs.

Given a geomval set, I need to make updates to each table. Here, I currently make use of st_setvalues (specifically variant 5), but I have to loop through each table to make those changes.

To speed this process (which can be quite slow, depending on the size of the geomval set), I'm thinking of making these individual tables as bands in one single table.

Does there exist a way to quickly update all bands at once so that I don't have to run the spatial intersection multiple times?

Currently, st_setvalues only operates on a one band at a time.


texts: Convert raster to vector using threshold values in an index - Landsat 8 SR

texts:
I'm looking to convert values between 4.45 and 10 in index FI into a vector polygon. I will be then be using this vector to clip my images prior to clustering. So the result must be clusters of L8 collection within the AOI defined by those said values in FI index. However, that is not the result produced. The vector conversion works but it includes two zones - i) between 4.45 and 10 and ii) outside these values.
Could someone help please? Code below.
var clusters = L8sr_ndvi_lswi_bsi_fi.map(function(image){

  var imgClip = image.clip(AOI)
  var imgFi = imgClip.select("fi");//
  //var zones = imgFi.updateMask(imgFi.lte(10)).updateMask(imgFi.gt(4.45))//
  var zones = imgFi.lte(10).add(imgFi.gt(4.45));
  zones = zones.updateMask(zones.neq(0));

  // Convert the zones to vectors.
  var vectors = zones.reduceToVectors({
    geometry: AOI,
    crs: imgFi.projection(),
    scale: 10,
    geometryType: 'polygon',
    eightConnected: false,

  });

  var updImg = imgClip.clip(vectors)//
  var training = updImg.select("B1" , "B2", "B3",  "B4" , "B5", "B6", "B7","lswi","ndvi","bsi")
  .sample({
  region: vectors,
  scale: 30,
  numPixels: 100
  })

  // Instantiate the clusterer and train it
  // specify the number of Clusters
  var clusterer = ee.Clusterer.wekaKMeans(10).train(training)
  return imgClip.cluster(clusterer);//result =
})


texts: Select by location using separate features iteratively (pyqgis)

texts:
I am looking for a way to run select by location, specific to the features of the 'INTERSECT' layer. I have got bus locations and bus stops and want to figure out when the bus has entered and left individual stops.

Here the orange dots are the bus stops (with varying radius) and the purple ones are the GPS locations of one bus.

With the standard logic, select by location select the GPS locations in all stops and I cannot discriminate when the bus has visited individual stops.

I am using this command for the selection:

processing.run("native:selectbylocation", {'INPUT':Course,'PREDICATE':[0,6],'INTERSECT':Stops,'METHOD':0})


texts: Convert shapefile > 4 GB to ESRI Geodatabase

texts:
I have 2 large shapefiles > 4 GB and 2 GB, is it possible using QGIS or ArcGIS to convert this to an ESRI Geodatabase so I can work with the data?

I have ArcGIS Advanced License 10.5.1 and QGIS 3.2.

I try adding them to ArcMap- right-click layer - Data - Export data, select file geodatabase as output but as far as I can see Arc doesn't support that file size for .shp files. So when I load them in, it doesn't load more than the 2 GB, which is also the only part that is converted to the Geodatabase. Hence I have a loss of data when converting that way.


texts: Join command, desired field not available to select

texts:
I am using ArcGIS Desktop 10.1.

I want to join a shapefile in which I stored geographical information with a CSV data file containing certain variables. However, when I try to apply the "join" command, the variable from the CSV file (6 digits zip code) to base the join on is not shown. Therefore I have exported the CSV file into a DBF file, but the variable still won't show up in the join command. For your information: The variable is shown as "Type:Long".

What is the solution to this problem?

I have read somewhere on the site to "Add a field" and then use the "Field Calculator" and simply copy the column in question with a Python command (!!). I tried this and copied the column, and now it shows "Type:String", but the join is still not functioning.


texts: Google Earth Engine: Converting raster to vector in the same resolution as the original image

texts:
I am using Google Earth Engine to reduce a raster region to polygon vectors. What I want to achieve is:

Obtain a set of polygons each of which matches the grid of measurement in the raster dataset. If the raster dataset's measurements are 0.25 x 0.25 degree, I want to break up the region in the map to polygons that align with these 0.25 x 0.25 degree grids.

Here's how I am doing the raster to vector conversion:
// a "geometry" polygon is imported

var raster_features = ee.ImageCollection('ECMWF/ERA5/DAILY')
                     .select('maximum_2m_air_temperature')
                     .filterBounds(geometry)
                     .filterDate('2017-01-01', '2018-12-31');

var mean_features = raster_features.mean().multiply(10000).toInt();
var scale_to_use = mean.projection().nominalScale();

var vector_features = mean_features.reduceToVectors({
  reducer : null,
  geometry : geometry,

  //scale: 28000,
  scale: scale_to_use,

  crs :mean.projection().getInfo().crs,
  //crsTransform :mean_features.projection().getInfo().transform,

  geometryType: 'polygon',
  labelProperty : 'grid',
  eightConnected: false
});
Map.addLayer(vector_features, {}, 'vector_features')

The resolution of ECMWF/ERA5/DAILY dataset is 0.25 arc degrees. I was expecting the length and breadth of the grids returned by reduceToVectors to be 0.25 deg, but they are 1 deg instead. I want to get polygons in the resolution that the dataset is in, which is 0.25 degree. How do I achieve that?


texts: Selection on Feature layer (geoprocessing)

texts:
I have a feature layers with states and counties fields. I want to make a selection tool. Tool dialog should have input for feature layer and two drop downboxes for states and counties. When user gives input feature layer, states drop downbox should get updated with unique values fron states column so that user can select state. And when user select state the counties dropdown box should get updated with the counties of the selected state, so user can select county from the selected state only. After that when user click on ok button the respective county should get selected on map.
Can anyone please tell me how can I achieve this using model builder or python script in arcgis desktop?


texts:
I am using Huawei MediaPad M3 Lite (Android 7.0). QField suits my requirement.
However at time I am unable to close the panel on the right after I have viewed the attributes of the layers.

The other problem is that I cannot snap pictures and position it on the map.

How can I rectify this problem because I need these to function well?


texts:
I'm trying to run QGIS and I am getting constant minidump crashes and "qgis-bin.exe has stopped working" crashes.

I have installed QGIS 2.18.0, the 64-bit version, on a 64-bit Windows 7.  The installation went smoothly except for the note that it "Could not find C:/Program Files (x86)/postinstall.log".  But when I try to add an ESRI shapefile or raster layer the program crashes with one or the other of these two errors.

This happens in QGIS Desktop both with and without the Grass 7.0.5 add-in.  Before installing 2.18.0, I had 2.16 installed for a couple of weeks and the same thing was happening.  I uninstalled 2.16 before installing 2.18 the first time.  I have uninstalled and reinstalled 2.18 to no avail.

I'm trying to add layers from the local drive, not an external source.  I've tried right-clicking, double-clicking, and dragging the layers into the layers panel and into the map area.  I've also tried the Add Vector Layer and Add Raster Layer buttons.  I've tried adding the sample data included with the download.  I've tried unplugging all the plug-ins first.  I cannot progress even one step before it crashes.

In contrast, the QGIS Browser module seems to work just fine.

Is this likely to be a problem with the install, or with the way I'm using the program?  I'm a newbie with QGIS so maybe I'm just doing something dumb.


texts: Incorrect creation of flow direction raster from interpolated DTM in ArcGIS

texts:
I have DTM for quite small area (400x200 meters) which was interpolated from elevation points. I want to create raster which will show water stream (not arrows but basically slope and river flow map - output from Flow Length would be enough). My DTM looks correct.
I use Fill to fill gaps and then Flow direction, where the problem occurs as you see below. For some reason there are created these straight lines which then also influence the Flow Length output. Please do someone know where is the problem and how to avoid it? Maybe the problem is that the DTM was created by interpolation, the area is to small or that there is small "lake" in the middle.


texts: What happened to Equidistant Conic Projection in Oracle 12c?

texts:
Back in 10x and 11x (I think), you could use Conic (sp, Simple / Equidistant) projections, but working in 12c, I can't figure out if they are still supported and if so, where to find them.

I see what Oracle was trying to do in making the CRS support more flexible, but to me, it is just incomprhensible.

At a basic level, I want to impletment EPSG/ESRI 102005 in Oracle.


texts: Selecting only visible layers exluding values that have not been symbolized?

texts:
I have many mxd's and I want to capture the features that are visible in a specific layer and feed that to another data processing tool.  I wrote some code that gets a data frame extents and selects the features in the given extent. This seems to work pretty well.  However the issue that I need help with is figuring out how to only select what is displayed.  For example, a feature class might be symbolized by category, and have an " symbol value that is not checked, see screen shot.

Is there a way to select only the features from symbols that are visible?


texts:
I have used Osmosis to load a large cut of OSM into a PostGIS DB, that works great. BUT, I would like to use the OSM-Bright style (carto CSS) to style it in TileMill and built tiles, offline. Going through the OSM-Bright user guide/install etc. there is no option to configure it to use an Osmosis schema (only osm2pgsql and imposm).

On windows osm2pgsql will not import very large OSM datasets hence the Osmosis route to get me this far.

Is there a way to style my data imported with Osmosis, or modify my DB so it is in the correct schema (the schema osm2pgsql uses)?


texts: Refreshing Catalog window from ArcMap using Python?

texts:
I am running my script from ArcMap and I have the Catalog window open on the side. I delete some extra files in the end of my script but I have to refresh the Catalog window to see them go. So I decided to use the following code to refresh it automatically:

import arcpy
from arcpy import env
import os.path

# Parameters
XS = arcpy.GetParameterAsText(0)

# Environment Definition
desc = arcpy.Describe(XS)
path = desc.path
env.workspace = path #Everything is saved at *path*

arcpy.RefreshCatalog(path)

But it doesn't work either. Any idea what is missing here?

I am running my script though a toolbox.
And my files are located at path.


texts: Save as .csv with coordinates in both QGIS 2.2 and PyqGIS

texts:
These questions have been asked before, but since they are not answered I will try to pick up the ends here.

1) First, as stated here:

QGIS 2.2 saving layer as CSV with GEOMETRY=AS_WKT not working

in Qgis 2.2 it is not possible to save a shapefile as .csv including the coordinates using the function _GEOMETRY=AS_WKT_. The issue

http://hub.qgis.org/issues/9597

is set to closed. Does it mean that it has been fixed? And if so, how to make it work again?

2) Second, if you want to save a shape file as .csv with coordinates through PyQgis it can be done with:

https://gis.stackexchange.com/questions/61857/how-to-write-vector-layer-to-csv-file-with-geometry-by-python

as stated here:

How to write vector layer to csv file with geometry by PYTHON?

But again, this does not save the coordinates, and I can't find any solution to this either. So right now, I still have to use an old version on QGIS.

IS there any solutions to these, available at the moment?


texts: How to extract GeoJSON vector layer coordinates from DOM, Openlayers 3

texts:
I created a GeoJSON vector layer using OpenLayers 3 API, the vector layer is a county map. My goal is by clicking, I will be able to select the county and extract the coordinates. (like shown below)

My current method of extracting the coordinates is once on.select, goes to the DOM and find the coordinates. Something like this:

var coords;
var select = new ol.interaction.Select();
    map.addInteraction(select);

// var selectedFeature = select.getFeatures();
selectclick = select;
selectclick.on('select', function(e)
{
  if (e!== null)
  {
  console.log(e),
  coords = e.selected[0].c.target.s;
  console.log(coords);
  }
});

This way works find, but I had to manually go the DOM and find the coordinates myself, then write the location in variable "coords". This probably not a good solution when face dynamic environments. Is there a better way to find coordinates in my case?


texts: I am not a MapInfo user, but I have been asked whether that software can plot geotagged photos and link the images to the points (eg, click on a point and the photo pops up)?


texts: Deleting duplicate values (depending on second condition) from raster attribute table?

texts: How can I automaticly delete rows from a raster attribute table if the value of a specific column appears twice or more? For example my attribute table (below) shows in the column "Streamlink" the value 22 twice, now I want to keep the row with the maximum value in column "FLA_FLD666" and delete the smaller one. The best would be a arcpy solution.