texts: Voronoi polygons that run out to infinity?

texts:
I have been working on implementing the voronoi.py script (Bill Simon's translation of Steve Fortune's code).  According to the QGIS link it is basically the same code that I am working with, but the way to generate the polygons that can extend to infinity differs between the two.  In the QGIS case, I presume that some bounding rectangle is generated at some extent beyond the site point extents.  In my case, I generate false site points as a circle surrounding the site points of interest.  The circle points have a radius twice the diameter of the extent of the site points.  The results between the two are shown in the attached figure, QGIS implementation on the left and my implementation on the right.
.

As one can see, the QGIS implementation produces weird polygons which have no site points (see site 473 as an example in both implementations).  This is obvious wrong, but I am not sure whether my implementation is correct either.  Since I don't have access to an ArcInfo license for ArcMap, I can't compare it what ArcMap generates natively and I can't see a way to a way to attach the site point shapefile for someone to test.  So if anyone has any commentary on my use of an "infinity circle" to help force edge polygons versus an "infinity rectangle", I would appreciate it.  Better still, if anyone is interested in testing the output on some other GIS platform I can provide the shapefile of the site points I used.

It appears part of the problem is that the polygons beyond the site points overlap, so when symbolized without a fill, then appear as distinct polygons.  The following figure shows what happens when the polygon about site 473 and 415 are selected.  So the original questions remain, but there is a new wrinkle.

It is a bug (or an "incompleteness"), the correct results in QGIS if a buffer of 0% is used.  What still needs to be resolved is whether a dense pattern of circular "infinity points" will yield the same result as rectangular "infinity edges".


texts: Using PyQt5 Objects inside QTabWidget and QScrollArea - QGIS Plugin

texts:
I'm using PyQt5 and QT Designer to create the user interface for my QGIS plugin. Now I'm learning to use QTabWidget and QScrollArea. What I don't understand is how to use the objects inside tabs or scroll areas. For example in my dialog there are two tabs, inside these there's a scroll area and finally the other objects like labels, lineEdits and buttons. You can see this in the image below. The problem is that if I try to use, as usual in my run method, simple functions like:  self.dlg.label_2.setText() or self.dlg.lineEdit.text(), I receive this kind of error log: AttributeError: 'Dialog' object has no attribute 'label_2'.
How to solve it?


texts: Is it possible to save selected objects such as map, legend, scale bar be saved as image without saving the whole layout in QGIS?


texts:
Can someone help me on how I can switch the basemap layers?

I'm using Mapbox layers, but when I switch from satellite to road, the feature layer is pushed to the behind the new basemap.


texts:
I have a question concerning the SRTM elevation data (http://srtm.usgs.gov/).

As a little thought experiment, assume I want to build a 3D model of the earth's surface as a point cloud. The earth's center is at position (0, 0, 0). For each latitude/longitude position for which elevation data is available, I want to find the Cartesian coordinates of the point on the earth's surface.

It says on the SRTM website "Heights are in meters referenced to the WGS84/EGM96 geoid". I'm not 100% sure if I understand this correctly.

From what I understand, the meters in the SRTM data are above sea level with respect to the EGM96 geoid. So for certain coordinates (latitude, longitude) I first need to find the Cartesian direction vector for these coordinates (let's call it d). Then I need to find the sea level at these coordinates with respect to the EGM96 ellipsoid (let's call it s) and then the height above sea level for these coordinates from the SRTM data (let's call it h). Then the position of the surface point would be (h + s) * d. Is this correct or did I misunderstand something?

If this is in fact correct, I have two more questions:

How can I find the Cartesian direction vector from the latitude/longitude?
How can I find the sea level with respect to the EGM96 geoid at a certain latitude/longitude coordinate?


texts: How to interpolate GPS Positions in PostGIS

texts:
I have a PostGIS table of GPS positions for every five seconds:

2011-01-01 00:00:05, POINT(x1,y1)
2011-01-01 00:00:10, POINT(x2,y2)
2011-01-01 00:00:15, POINT(x3,y3)
...

I'm looking for a query that will return values (timestamp and point) for every second. It's ok to assume that points are connected by a straight line.

I'm specifically looking for a way to do this inside the database and not by writing some external script.


texts: QGIS show labels in beggining or ending points of Linestrings

texts:
I've got a Linestrings layer.

Each Linestring has among other fields, these ones:

  id  ENTRY EXIT
    1   S1  E2
    2   S1  E3
    3   S1  E3
    4   S1  E3
    5   S1  E3
    6   S1  E1
    7   S1  E3
    8   S1  E2
    9   S2  E5
    10  S2  E3
    11  S1  E1
    12  S1  E1
    13  S3  E4
    14  S1  E2
    15  S1  E1

I want to show for each Linestring, the ENTRY value as a label at the beginning of the Linestring and the EXIT value at the ending.

Something like this (Only the red labels):

One approach would be to create a new points layer with the beggining and ending points of each Linestring and add the data of the Linestring to each point and after that just show the labels...but, is there a way of achieving what I want using only the existing Linestring layer?


texts: Singlepart to multipart not producing results

texts: I have a line - a street downloaded from OSM.  It is in separate parts geometrically, but is contained in one record in the attribute table.  That means it is only one row.  I would like to split this street into its separate parts (you can see that it is made up of 4 separate lines), but I get this error when I try to run it through Vector > Geometry tools > Singlepart to Multipart.  Why is this happening? Why do I need more than one feature to run this tool?


texts: Using time to calculate service areas in ArcMap?

texts:
I am having trouble computing some service area polygons in ArcMap 10.2. It seems to work fine when calculating distance, but has problems when looking at time.

I have UK road data (Meridan 2), and am attempting to create drivetime polygons around a number of retail sites. For each road length Ihave assigned a speed (based on the type of road - ie major/minor and the settlement type - ie rural/suburban/urban) and calculated the time in seconds to travel along that length of the road.
I have created polygons using 'make service area' in the network analyst tool. I first looking only at distance (metres), and this worked fine (see image). However, when I try to calculate service areas based on time (seconds), I get very odd results.

I have copied below the service areas that I get around 5 sites, with an impedance of either 5000 metres or 0.2 seconds (obviously I want a longer drivetime than that, but anything larger than this gives huge service areas - this is the problem).

As you can see, when calcuatingareas based on time, some of the sites have tiny service areas (which makes sense, given the small time available) but some (one in this example) have very large areas. Actually in this example two of the service areas are even rectangles, which is even more surprising and I hadn't noticed before.

All street lengths in my network seem to have positive integer values in the 'seconds' field. These values for seconds tend to be quite similar to the values for meters that are working fine (all my speeds are between 6 and 30 metres per second). Seconds is specified as a cost attribute in the network dataset properties for my network dataset (all_roads_speed_ND). I really don't know what's going on.

Service areas with distance

Service areas with time


texts: Using "Refactor Fields" to Change Precision

texts: I am trying to use the Refactor Fields tool to change the precision of one of the Real fields in a shapefile from 0 to 4. I simply change the precision from a 0 to a 4 under the Fields mapping in the tool. After the tool runs, the field has been changed from a precision of 0 to a precision of 5, not 4. Not only that, it has changed the precision of every real number to a precision of 5 even those which had a precision greater than 5. Even if I don't change anything but just run the tool, it changes every real number to a precision of 5. Is there something more I need to do?


texts: How to add raster data to Table of Contents using python?

texts:
I can easily add shapefiles to the Table of Contents in Arcmap using the following codes:

mxd = arcpy.mapping.MapDocument("current")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
path = os.path.abspath("shapefile_path")
addLayer = arcpy.mapping.Layer(path)
arcpy.mapping.AddLayer(df, addLayer)
arcpy.RefreshActiveView()
arcpy.RefreshTOC()
del mxd, addLayer

but I am confused how to add raster data (which usually consists of several files). Any help would be much appreciated.


texts: Install Postgresql / Postgis with SFCGAL

texts:
I need to access ST_StraightSkeleton() function. At the moment i have a Postgresql
instalation PostgreSQL 9.3.4 compiled by Visual C++ build 1600, 64-bit installed via EDB with Postgis 2.1.7 on a windows machine that seems not have SFCGAL installed.

  1st What would be my chances of having SFCGAL installed on this installation ?

  2nd I'm willing to make a fresh install of another PostgreSQL / Postgis installation  (built-in with SFGCAL) if that is less effort, but im a bit confuse on which installations and the order of installation for the .exe files at Postgis site below

http://winnie.postgis.net/download/windows/pg93/buildbot/

EDIT 1:


texts: ArcGIS 10 using ArcPy to save feature class table as a geodatabase table (i.e. without geometry)

texts: I have a point feature class where I need to - in arcpy - save the feature attribute table i.e. without the geometry. I need arcpy to do the same as Opening Attribute Table then Export the table, but not with cursors as this is a very large attribute table (very many points).


texts: What is the maximum theoretical refresh rate for GPS

texts:
I understand that most devices have 1 Hz refresh rate for GPS positioning. I tried to read about the actual GPS implementation but I could not figure out how often the location can be determined.

Given a custom chip, how often (Hz) the real position of the GPS receiver can be determined using only GPS signal? GPS signal chip (/bit) length is around 1 µs. Would it be possible to compute location every 2-3 chip (/bit) if given enough CPU power (that is, compute again for each time signal changes between 0 and 1)?


texts:
I'm trying to style a torque heatmap in Cartodb based on a weighting value(pntcnt) rather than sum of points in an area. I think i've generated the appropriate value with the -torque-aggregation-function:"sum(pntcnt)";

However I can't seem to get the styling right.

I am new to CartoDB.

Full CSS code.

Map {
-torque-aggregation-function:"sum(pntcnt)";
-torque-resolution:2;
-torque-data-aggregation:linear;
}

#heatmap{
  image-filters: colorize-alpha(blue, cyan, #008000, yellow , orange, red);
  marker-file: url(http://s3.amazonaws.com/com.cartodb.assets.static/alphamarker.png);
[value > 500] { marker-fill: blue; }
[value > 1000] { marker-fill: cyan; }
[value > 1500 an] { marker-fill: yellow; }
[value > 2500 { marker-fill: red; }
}

Data look something like this

Cartodb_ID        the_Geom                  pntcnt
7384             -107.6227, 24.4173         3391
7476             -107.6140, 24.4219         2514
6481             -107.7091, 24.4809         2352
6572             -107.7004, 24.4764         2258


texts: Intersect Tool - Same Polygon Feature Class Only Want Inner Overlap

texts: I'm trying to find the best way to only get the inner overlap of multiple polygons within the same feature class.  I've tried the Union/Dissolve method as well as the Intersect tool and can't get the desired result I want.  My sense is that I may need to create some kind of topology rule to do it, but ostensibly I just want to find the inner areas where the polygons overlap and right now the software is selecting ANY polygon that overlaps with another. I've attached an annotated picture to show my challenge.


texts: Update table with condition from another table in postgresql

texts:
I have two table name servers and sites. Servers contains timing and workmode and geometry values. Where sites contains geometry values and an empty time_stamp column. I want to fill this timing column.
For this, I want to check for all servers data whether it is within 9 unit distant for each sites. Hence, servers data 1 will search for all the sites and if it finds a site within 9 unit then it checks for its own work mode,
if it is true then it writes the time in time_stamp column of that sites record. All the servers data with true work mode should find out one site. But there can be more sites. So, all the sites records might not update. The algorithm is like this:

FOR each server data
    FOR each sites
        IF current server data is within 9 m distant of current site AND workmode is true THEN
            sites.time_stamp <- server.timng
            break sites loop

For this I have written the following query. But not sure whether it is serving my purpose properly or not. Can anyone suggest? Should I need to add any loop here in my query?

    UPDATE sites st
SET time_stamp = timing FROM servers sv
WHERE st_distance(st.geom, sv.geom) < 9 AND sv.workmode = TRUE
;


texts: mapnik error - fatal flex scanner internal error--end of buffer missed

texts:
I am using Mapnik 2.0.1rc0 on a Windows 7 with python 2.7 32bits inside a GeoDjango App. I often got this error: fatal flex scanner internal error--end of buffer missed that makes python to crash (sometimes python crash without error message) when zooming in and out of mapnik tiles rendered in OpenLayers map. The tiled data are stored in spatialite 4.0.

I know such setup could provide a lot of source for this crash. Anybody has encountered this error before or may know where to look to make it works without crashing?

Edit:

GeoDjango is used as a Tile Map Server. I start by creating the xml in python then I use mapnik python bindings to generate the tiles. The data is already stored in spatialite. So each time I zoom in and out in the map, mapnik query spatialite like this:

    query = '(select ' + geometryField + ', "id" from "shapefile_feature" where' + ' shapefile_id=' + str(shapefile.id) + ') as geom'
    datasource = mapnik.SQLite(file='C:\mygeosite\sqlite3\sql3.db',
                                    table=query,
                                    srid=4326,
                                    geometry_field=geometryField)


texts:
I have written a code to check if areas will be inundates when a certain value is added, the code doesn't give any errors but at the same time it doesn't do anything.

Just says running code complete, I am not even sure what the problem is.

My code is here:

import arcpy
import os, sys, traceback, shutil
import numpy as np

arcpy.env.overwriteOutput = True

try:
     outputfolder = arcpy.GetParameterAsText(0)
     fc = arcpy.GetParameterAsText(1)
     inputraster = arcpy.GetParameterAsText(2)

     #check to see if the arcinfo license is on
     if arcpy.CheckExtension("Spatial") == "Available":
          arcpy.CheckOutExtension("Spatial")
     else:
          arcpy.AddError("Spatial analyst extension is not licensed.")
          sys.exit()

     #set the Progressor
     #arcpy.SetProgessor("step","Checking input requirements....",0,9,1)

     #Process :create an out put folder
     arcpy.CreateFileGDB_management(outputfolder,"output.gdb")
     outputgdb = outputfolder + "\\output.gdb"+ os.step

     #create path to store raster data
     clippedraster = outputfolder + "clippedraster"
     flooded = outputfolder +"floodedarea"

     #check to see if the inputs are in the same projection
     desccs_ras = arcpy.Describe(inputraster).SpatialReference
     desccs_fc = arcpy.Describe(fc).SpatialReference
     if desccs_ras.name == "Unknown" and desccs_fc.name == "Unknown":
          arcpy.AddMessage("Your data does not have a defined spatial reference")
     else:
          if desccs_ras.projectionName == desccs_fc.projectionName:
               ""
          else:
               arcpy.AddError("The input Data should be in the same projection")
               sys.exit()
          del desccs_fc,desccs_ras

     #remove the areas below sea level
     currentRaster = arcpy.Raster(inputraster)
     rastertoclip = Con(currentRaster>=0,currentRaster,1)

     #clip the raster data to the interested zone extend
     arcpy.Clip_management(rastertoclip,"#",clippedraster,fc,"#","ClippingGeometry","MAINTAIN_EXTENT")
     clip_raster = arcpy.Raster("in_memory")
     areasbelow = clip_raster <= 2
     areasbelow.save(flooded)

     md = arcpy.mapping.MapDocument("CURRENT")
     df = arcpy.mapping.ListDataFrames(md)[0]
     result = arcpy.MakeRasterLayer_management(flooded,"FloodedArea")
     layer = result.getOutput(0)
     arcpy.mapping.AddLayer(df, layer,"AUTO_ARRANGE")

except:
     tb = sys.exc_info()[2]
     tbinfo = traceback.format_tb(tb)[0]
     pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n    " + \
            str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"


texts: lulc mapping from modis surface reflectance

texts: can i prepare LULC with 10 classes from MODIS surface reflectance level 3 data 500m resolution... if so what are the preprocessing techniques should i do on the data..i have to prepare LULC for yearly 2000-2014 for western ghats on a regional scale . will the size of training sets alter the accuracy of LULC? i tried Maximum likelihood classification with green,NIR,red combination and there was a lot of variation from year to year and the map was not accurate. i didn't process the data as it was already processed level 3 .


texts: Connecting PCI, ArcGIS, and XTools with Python?

texts: I have a process that I'm trying to fully automate. There are 3 broadly defined chunks. An intial ArcGIS chunk that uses XTools. A PCI chunk. and a plain ArcGIS chunk. PCI python requires 64bit arcpy in order to connect. Is there any way to do this and still use XTools, which I understand, require 32bit arcpy?


texts: Making a map with 8 different rasters of the same place but with only one legend in QGIS

texts: I have 8 raster images of a farm that are changing with respect to time. I want to create one map as shown in the attached image which will show the 8 raster images with only one legend to represent raster pixel values across all 8 rasters.


texts: QGIS doesn't create the correct projections when creating new shp files

texts:
I received shapefiles exported from ArcGIS with EPSG:2227 as the CRS, which should be:

  +proj=lcc +lat_1=38.43333333333333 +lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000.0001016 +y_0=500000.0001016001 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

When I loaded it into QGIS, it says that it is a user defined CRS with the lat_1 and lat_2 mixed up and they have slightly different x_0 and y_0 values:

  +proj=lcc +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000 +y_0=500000.0000000001 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

Even when I reproject the .shp from the user CRS to EPSG:2227 and then load it back into QGIS, it says that it is a user-defined CRS.

Has anyone else seen this happen?


texts: how to get the workspace names from geoserver

texts:
i am a student and i am trying to build a web application based on geoserver. I have geoserver installed on my PC and i have already build an html multi-form.

I want when user selects option A, this option to be linked with my A geoserver's workspace, and when user selects option B, to be linked with B workspace.

After the selection a new form appears (i have made it with jquery) and i want whatever selection be made above (main form), my new form to get the layerlist from the selected workspace.


texts: Area of the city, I have amenities, roads and and admin relations

texts:
I am a complete beginner in spatial data and I do not know how to calculate the area of the city having amenities, roads and admin relations in postgres database. I tried searching the internet, but the answers I found are based either on the contour of the city, or a given polygon, or postal codes. I do not see a possibility of using these solutions in my case.
Here is my data, the city I am looking for is called Kraków (of Cracow). Here is my data in ESRI shapefiles format.

https://drive.google.com/file/d/12snSgNc4SrtG6x0HeBMFIzvYH8Wfsl_j/view?usp=sharing

Here is the same data but with created sql files.

https://drive.google.com/file/d/1opyg0njPHI2L2jlrlQ6nnQmxrkCijYao/view?usp=sharing

I will be very thankful for helping me what i should do. I bought the courses on Udemy, but it is going to take some time before I watch all of them. I know that I should use ST_Area function, but I do not know what to put as arguments. I should probably have a polygon, but how to create it? I have points, but how can I know which point belongs to Krakow, if there is no column postal code.


texts: Calculating Flow Direction using ArcMap with mixed X,Y,Z units

texts:
I have a DEM in UTM that I am converting to a StatePlane coordinate system in feet. This raster is going to be used to generate watershed delineations - therefore I will generate a flow direction raster from the DEM data. I don't think it matters (according to my understanding and some back-of-the-envelope) examples, that the z-units are in meters and the linear units are in feet, but could it in some edge case I can't think of?

How do you generally handle operations on a DEM with mixed units? Do you always convert the Z-units?


texts: Creating a layer from a selection, when I try to symbolize or label the layer, it brings back all the non-selected features

texts:
I'm trying to create a layer from a selection in another layer in ArcMap 10.5.1 like so...

Doing Right Click Layer - > Selection - > Create Feature from Selected Features after doing my Select by Attributes and Locations. Those are just a Select by Attribute of this PREFIXDIR = 'W' AND STREETNAME = 'WILLOW' AND STREETSUF = 'ST' and then select by location all locations like such...

Then creating the new layer and showing that...

All is great! But as soon as I try to label the features/symbolize like the previous layer...it brings back all the non-selected features.

Using the following Python expression to label...

And the symbolology I'm just importing from the previous layer.

As soon as I label or import symbology all the previous features are brought back...

Like for example when I import the symbology

Anyone know what I'm doing wrong?


texts: Grouping wildfires dates by season and getting sum and count using ArcPy?

texts:
I am having some trouble creating a new shapefile and CSV to group wildfires into seasons.

What I want to do is to group the wildfires by meteorological season (summer, fall, winter and spring) for each year. Then, I would like to know how many wildfires and the area burned in each season for each year from 1950 to 2017.

Summer starts from 1st of June to 31st August, fall from 1st September to 30th November, winter from 1st December to 28/29 February and spring from 1st March to 31 May.

The result might be similar to this by creating a new field called "Season", another one called "Count" (to see how many wildfires are in that season) and another column called Season_acres which will be the sum of all acres burned in a season.

So far I have got these two separate scripts:

Using an update cursor

    fc = "fires_1950-2017.shp"
    newfield = "MET_SEAS"
    fieldtype = "TEXT"
    fieldname = arcpy.ValidateFieldName(newfield)
    print "cleaned up fieldname:", fieldname
    arcpy.AddField_management(fc, fieldname, fieldtype, "", "",12)
    print "New field created!"
    delimfield = arcpy.AddFieldDelimiters (fc, "ALARM_DATE")
    cursor = arcpy.da.UpdateCursor(fc, ["MET_SEAS"], """ALARM_DATE >= date '1950-06-01 00:00:00' AND ALARM_DATE <= date '1950-08-31 00:00:00'""")

    for row in cursor:
        row[0] = "Summer 1950"
        cursor.updateRow(row)
    del cursor

Using Statistics analysis, which I am not able to use with an update cursor

    in_table = "fires_1950-2017.shp"
    out_table = r"D:\edited_data\try.csv"
    stat_fields = [["GIS_ACRES", "SUM"], ["GIS_ACRES", "MEAN"]]
    stats = arcpy.Statistics_analysis(in_table, out_table, stat_fields)
    print "script finished"


texts:
When I'm working on my routing, I wanted to make a topology and this error occur
saying:

  Unexpected error parse error - invalid geometry

this is my code

select pgr_createTopology('table_name',0.001,'the_geom','gid');

my RSID is - espg 5234

this is my geometry column with table


texts: Loading *.rvc files from TNTmips into QGIS?

texts:
I have many files obtained for a class I took in 2004. We used TNTmips to work with these files.

The file extensions are ".rvc". They show up as executable files in Mac Finder. I remember seeing somewhere that the TNTmips software used the ".rvc" file format for vector data, raster data, and -- something else.

It has been so long that I do not remember what our class used these files for, or how we used them.

Can files with the ".rvc" extension be loaded into QGIS? If so, how?


texts: ArcGIS Feature Layer from Feature Collection using ArcGIS API for JavaScript?

texts:
I am trying to create a polygon feature layer from a feature collection using the ArcGIS API for Javascript.  I have modified Esri's Flickr sample to use polygon geometries, rather than point geometries. The features are not getting added to the layer properly. I have put together a JSFiddle here that illustrates the problem. I have commented out the method where I add the features (line 137).

http://jsfiddle.net/gregKnight66/mmhqqan6/

There isn't an error per se, I am just not getting results (the features are not added to the layer).

var featureLayer = new FeatureLayer(featureCollection, {
    id: 'buildingFeatures'
});

map.on("layers-add-result", function (results) {
    addBuildingFeatures();
});

map.addLayers([featureLayer]);

function addBuildingFeatures() {
    var features = buildingFeatures;
    featureLayer.applyEdits(features, null, null);
}


texts: How to Mosaic Images using ArcMap when coordinate system is undefined?

texts: I have images which are not Georeferance with any coordinate system and i want to mosaic those images in ArcMap 9.3


texts: Is there any way to make QGIS layout to activate/deactivate certain layers?

texts:
I have a raster collection over certain area. Every raster is on top of another one. So if I have all of them visible, I can only see the one on top.
What I need is some kind of plugin (or feature I am not aware of) so for each layout in a workspace. I can select what layer are visible or not, and then somehow use a plugin like automapper to extract everything as PNG without having to go manually all by myself.
I don't mind coding a solution if it is not a thing out of there. In that case, I would like also some suggestion on where would be a good place to start.


texts: How to save a clipped NetCDF data as a new NetCDF file in Python?

texts:
I have used xarray to clip my NetCDF file (it can be downloaded here) by a shapefile (it can be downloaded from here). Here is the code:
import xarray
import rioxarray
import geopandas
from shapely.geometry import mapping

xds = xarray.open_dataset(
    "3B-DAY.MS.MRG.3IMERG.20170216-S000000-E235959.V06.nc4",)

xds = xds[['precipitationCal', 'precipitationCal_cnt']].transpose('time', 'lat', 'lon')

xds.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True)
xds.rio.write_crs("EPSG:4326", inplace=True)
geodf = geopandas.read_file("NAME_LONG_201.shp")

clipped = xds.rio.clip(geodf.geometry.apply(mapping), geodf.crs)

Now, I have a clipped object which is of type: xarray.Dataset. How can I save this dataset into a new NetCDF file?


texts: Creating model that simulates water flowing down crack, then splitting at 'Y' joint?

texts:
Maybe this is a simple problem, but I've been given a task of simulating how water flows in one line, along a crack, and then hits a 'Y' junction, where half the water flows in each of the split directions.  I'm having trouble imagining how to do this.

I had thought about just creating 2 routes, that diverge at the 'Y' joint.  Is there a better way to do this?


texts:
I have two shapefiles which should be connected to one single shapefile without a gap, unfortunately there is a gap between them:

I tried to use v.cleanin Grass GIS but because of the shape of the two lines, it was not possible to connect only the wanted area:

Does someone has an idea how to solve this problem? I also tried SpatiaLite in QGIS which also didn't work well.


texts: Is it possible editing data inside QGIS Server?

texts:
I need to modify some field, for example i have a polygon but i would like to add a picture or pdf to it, so if someone click on the polygon he could see the picture or a document.

Is it possible modify the data on QGIS Server?


texts: Why does QGIS save EPSG:4326 as long, lat?

texts: Why does QGIS save EPSG:4326 with coordinates in the longitude, latitude format? Isn't the correct standardized format latitude, longitude?


texts: Create a dbf table by copying the structure of another dbf table

texts: I have a dbf table and I wish to automatically create another dbf table with a different name but with exactly the same structure i.e.fields as the existing table. Is there any ready code in VBA and ArcObjects for that? Thanks


texts: Creating animated graphics to symbolize selections in ArcGIS Explorer?

texts:
I want to create and add some animated graphics to the features based on a query in ArcGIS Explorer, like some moving waves at the queried features.

Is it possible?


texts:
I am using a script written by John to import a csv file and make a point feature class. My table looks like this:

"timestamp_pretty","mmsi","imo","shipType","lat","long","month","week"
"31/12/1969 23:59:59",265650950,9544322,"HSC",57.661503,11.84175,12,1
"31/12/1969 23:59:59",305706000,9583902,"CARGO",59.486267,20.078266,12,1

Tthere is a date field in my table. How to import it as date format instead of text?

# Author: John K. Tran
# Contact: jtran20@masonlive.gmu.edu

import arcpy
import os
import csv
import time
import locale
from arcpy import env

start_time = time.time()

##incsv = arcpy.GetParameterAsText(0)
##outfc = arcpy.GetParameterAsText(1)

incsv = r"E:\DensityMaps\TestFilesLoop\test.csv" # Change this to the path of your CSV file.
outfc = r"E:\DensityMaps\DensityMapsTest1.gdb\TestSampleFC" # Change this to the path of your output FC.

spatialref = arcpy.SpatialReference(4326) # Create the spatial reference object as WGS84. Can modify if desired.

if not arcpy.Exists(outfc): # Create the output feature class if needed.
    arcpy.CreateFeatureclass_management(os.path.dirname(outfc), os.path.basename(outfc), "POINT", None, None, None, spatialref)

csv.register_dialect("xls", delimiter=",", lineterminator="\n") # Register the dialect for native CSV syntax in Microsoft Excel.
f = open(incsv, "r")
reader = csv.reader(f, dialect = "xls")

headers = reader.next() # Read the first line as the header names.
headerDate = headers.index("timestamp_pretty")

for header in headers[0:]: # Add fields for remaining columns if needed. Default is TEXT field.
    arcpy.AddField_management(outfc, header , "TEXT")

    print "Header: "+header

cursor = arcpy.da.InsertCursor(outfc, ['SHAPE@XY'] + headers[0:]) # Create InsertCursor.
count = 0
for row in reader:
    if count % 1000 == 0:
        print "processing row {0}".format(count)

    date = row[0]
    Ycoord = row[4] # Make sure 'Lat' is in the 4th column.
    Xcoord = row[5] # Make sure 'Lon' is in the 5th column.
    newrow = [(float(Xcoord), float(Ycoord))] + row[0:]
    cursor.insertRow(newrow) # Insert point in FC for each row in CSV.
    count += 1

del cursor
f.close()
print ("--- %s seconds ---" % (time.time() - start_time))


texts: Getting PNG map from PostGIS to Node and then jquery

texts:
Using Jquery, I am aiming to get png map from postGIS which has been created dynamically from node server in different port. I used jsonp to overcome cross domain issue in jquery.

$.ajax( 'http://localhost:3000/terrain', {
              type:'GET',
                  crossOrigin: true,
                  data:koordObj,
                  contentType: "content-Type: application/image",
                  dataType: "jsonp",
                  success:function (data) {
                  console.log(data);
              }
            })
        .error( function() {
            console.log("error");
        });

server creates the image as :
"
localhost:3000/terrain/?callback=photos&minx=4000660.554526325&miny=4862598.159515699&maxx=4014342.782589372&maxy=4869286.399490653&terrain=slp&zoom=14&_=1459928691680
 "

But firefox arises this error:
SyntaxError: illegal character
�PNG
What might the problem be?


texts: Fastest strategy for proximity searches in SQL Server

texts:
I am using SQL Server 2012.
I'm implementing a back end for a mobile app which will have to do proximity searches to find nearby POIs (points os interest). I know it's a very common scenario and looks very simple, but there are many different ways I could implement it, so I would love to see how more experienced professionals are implementing these simple spatial searches.
Since a POI is just a POINT we don't need any complex calculations involving intersections or the like. That's why I initially thought that using GEOGRAPHY columns and spatial indexes could be overkill or even slower than other strategies. So I've narrowed it down to 3 approaches:
1) GEOGRAPHY column + Spatial Index
This is perhaps the de facto solution to this problem. Since we have spatial indexes and geography columns we could just use it and search by distance. Something like this.
SELECT * FROM POIs WHERE Loc.STDistance(@radius) <= @distance;

Since we have a spatial index on Loc, it should be very fast.
2) Using a "bounding box" over Latitude and Longitude columns
This is the trivial approach without involving spatial indexes. We find a bounding box for our point and radius then simply search on the Latitude and Longitude columns. If both are indexed this search should be very fast. We'll have to apply the distance function to filter out some values outside the "circle" but withing the bounding box. But that should be pretty fast. This idea is better explained here: http://www.movable-type.co.uk/scripts/latlong-db.html
Something like this:
DECLARE @lat float
DECLARE @lon float
SET @lat = -23.001029
SET @lon = -43.328422
DECLARE @maxLat float, @minLat float, @maxlon float, @minLon float
DECLARE @R float
DECLARE @distance FLOAT = 100 -- A distance in meters
SET @R = 6378137 -- Earth
SET @maxLat = @lat + DEGREES(@distance/@R)
SET @minLat = @lat - DEGREES(@distance/@R)

SET @maxLon = @lon + DEGREES((@distance/@R/COS(RADIANS(@lat))))
SET @minLon = @lon - DEGREES((@distance/@R/COS(RADIANS(@lat))))

SELECT * from POIs
WHERE
        Lat Between @minLat And @maxLat
    And Lng Between @minLon And @maxLon

3) Use an integral GEOHASH stored on an indexed column
This approach is very interesting and it is something people are using together with REDIS ordered sets to do proximity searches. The principle can be transposed to SQL Server by using an indexed column that stores the integral GEOHASH.
I've got this idea from Ardb: https://github.com/yinqiwen/ardb/wiki/Spatial-Index
It's also explained in a little friendlier manner here: Using geohash for proximity searches?
In other words one would compute a GEOHASH with a bit-depth corresponding to the radius of the search one desires, then compute 8 neighbors geohashes and finally submit a search using these geohashes as bounding boxes on the indexed column. This will be 9 BETWEEN operators on the WHERE clause of the SQL... The results will have to be filtered out due to some spurious POI being returned.
But it looks that this will be slower than method 2 as the where clause will be more complex although it will only query over a single column instead of two.
Is there a better/correct approach to this?


texts: Unable to add CSV files as a hosted feature layer in ArcGIS online

texts:
I have a hosted feature layers of 2 geometry types.

Geometry Type - Point
Geometry Type - Polyline

I want to export these layers and import it on another map.
I used the following steps:

Exported these hosted feature layers as a CSV file.

Opened map and click on "Add" button and selected Add layer from file option then I selected downloaded CSV file.

I am getting the following error:

For Point type

For Polyline type

I have uploaded the downloaded CSV file on drive
https://drive.google.com/drive/folders/1OY5uxHsvuL1JnzmaeL5p5Yv6fgIy7fsI?usp=sharing


texts: Handling MapQuest different zoom levels for lower 48 USA and rest of the world

texts:
The MapQuest open data site states:

  Please also note that global coverage is provided at zoom levels
  0-11. Zoom Levels 12+ are provided only in the United States (lower
  48)

I want to be able to use the MapQuest 12+ zoom levels in the USA lower 48 and provide my own tiles for outside the USA lower 48.

I am using OpenLayers and with GeoServer for my own tilesets for outside the USA

What is the best way to determine if my viewport or map center is contained in or out of the lower 48 so I can switch layers?


texts: Making raster output based on comparison between 2 input rasters using ArcGIS for Desktop?

texts:
I have 2 input rasters of rainfall for summer/winter and I want to make a 3rd output based on which one has the max value for that cell. So if the max value at the cell is summer it will assign 1, or winter 2 for the new output raster. That way I will have a raster showing which season has the max amount of rainfall across the different areas of the region.

For example if,

Summer(1) was,

200 300

300 350

and Winter(2) was,

300 350

350 200

The output would be,

2 2

2 1

I understand that the logic should be:
If max = summer
  output = 1
else max = winter
  output = 2

But I am unsure how to actually translate this into something like raster calculators conditional statements.

I am currently using ArcMap 10.2


texts: Creating raster basemap for web application?

texts:
I'm sure this is probably a common problem and many have already found the right answer.

I have a large mosaic raster file (a city's ortophoto) in .ecw fromat  and we want to use it as a basemap in a web application. Due to its size (2GB) rendering is very slow, therefore increasing performance is vital.

1) My first option was importing the raster in our ArcSDE geodatabase (on postgreSQL)-> build pyramids-> then publishing it as a map service with a 5 scale cache on our ArcGIS Server. Failed for unknown reasons (so far).

2) The second option would be importing the raster in the ArcSDE GDB then->create pyramids-> publish it as an image service to ArcGIS Server.

3) Third option would be to keep the raster in the folder-->copy the raster to the server --> publish it as an image service.

4) Fourth option would be to save the raster as a layer file-> then publish it in a map service.

Which of these four options would be the best when talking about rendering speed?

I use ArcGIS 10.1 + POstgreSQL DB + ArcGIS Server Enteprise 10.1

As brenth was saying you cannot publish an .ecw. I do have all the tiff files that created the .ecw mosaic, so I will try to create a mosaic from them and then publish using the fourth method. Does it sound ok?


texts: Can I batch import attribute tables into existing excel spreadsheets

texts:
I am using arcGIS 10.2 and I need to import the attribute tables of 100 layers into existing excel spreadsheets which will make graphs based on the values in the attribute tables.

Doing this manually is as simple as copying and pasting the attribute tables into the excel spreadsheets. However, this is quite tedious.

Is there a way that I can automatically insert these tables into the right sections of the excel spreadsheets?


texts:
I am trying to create a python scrip to batch erase shapefiles.  I have created a script that works, however, I know there is a better way and I just am not good enough at scripting to figure it out.  I want to create a loop to allow the script to just to go through the folders and just erase the new shapefile from the old to see change.  The way I have it now I have to update about 16 path names.  Also, the the script will send me an error if there is nothing to erase.  I know this requires an if statement but again, I am just now knowledgeable enough to figure it out.  Below, is the scrip that works and hopefully it will give and idea of what it is I am trying to achieve.

import arcpy
#Set workspace
arcpy.env.workspace = "USER_TEMP\\RAY\\test\\"

#Erase workspace
erase_ADE01 = "\\ADE01_ALL_BUFFERS_05_09_16.shp"
erase_COF01 = "\\COF01_ALL_BUFFERS_05_09_16.shp"
erase_CVI04 = "\\CVI04_ALL_BUFFERS_05_09_16.shp"

#Current buffers workspace
current_ADE01 = "\\ADE01_ALL_BUFFERS_TEST.shp"
current_COF01 = "\\COF01_ALL_BUFFERS_TEST.shp"
current_CVI04 = "\\CVI04_ALL_BUFFERS_TEST.shp"

#Outputworkspace
output = "\\test\\delta_test\\"

#Erase tool; NEED HELP-WHEN THERE IS NOTHING TO ERASE SCRIPT FAILS
arcpy.ET_GPErase(current_ADE01, erase_ADE01, output + "ADE01_ALL_DELTAS_TEST.shp", 0.0)
arcpy.ET_GPErase(current_COF01, erase_COF01, output + "COF01_ALL_DELTAS_TEST.shp", 0.0)
arcpy.ET_GPErase(current_CVI04, erase_CVI04, output + "CVI04_ALL_DELTAS_TEST.shp", 0.0)

print "Script completed"

UPDATE:

So, I don't know why I didn't think of this before but what I did was merge the erase features together so that way there is only one feature to erase from.  Then I just threw it in a "for loop" and everything works well for that part.  Now I am still getting that error for when there is nothing to erase (the features are the same).  How would I write an if statement for this:

#Erase tool; NEED HEELP-WHEN THERE IS NOTHING TO ERASE SCRIPT FAILS
for fc in arcpy.ListFeatureClasses():
   erase_feature =  "\\all_erase.shp" #this is the merged feature
   erase_output = "\\test\\delta_test\\" + "_ALL_DELTAS_TEST.shp"
   arcpy.ET_GPErase(fc, erase_feature, erase_output, 0.0)

UPDATE2:

So, I am getting closer to figuring this out.  My script will run all the way to the end now.  However, it wasn't erasing because of memory issues.  So i upgraded to python64x.  Here is my script:

for fc in arcpy.ListFeatureClasses():
prefix = fc[:6]
erase_output = "\\test\\delta_test\\" + prefix + "_ALL_DELTAS_TEST.shp"
try:
    arcpy.ET_GPErase(fc, merge_output, erase_output, 0.0)
except arcpy.ExecuteError:
    print arcpy.GetMessages(2)
    print "Erase for " + prefix + " unsuccessful"
except:
    print "There has been a non-tool error"

The "There has been a non-tool error" print pops up now that I have downloaded the different pywin.  Before I did that, the error was in the first except because of the memory issues.  What am I doing wrong?


texts: Moving polygon feature to be centered on given point using PyQGIS

texts:
I am trying to automate a process with inputs of:

Multipart polygon geometry
Four points - which are the actual center points of the four small rectangles within the larger polygon

How can I move the polygon to align with the four points using PyQGIS?
My first attempt was to start with the center points and construct the polygons around them, calculating the orientation of the points and rotating the polygons accordingly (see below), but I'll be receiving different polygons (with varying complexity) and hard-coding the dimensions will be too cumbersome. My next thought is to calculate centroids of those inner rectangles, but then I don't know how to snap them to the actual center points. I don't necessarily need the code, but a general approach would be helpful.
layer=my_layer
if not layer:
    print ("Layer failed to load!")
#loop through features
points = layer.getFeatures()
metrics={}
d=QgsDistanceArea()
for point in points:
    #calculate pairwise distances and bearings
    compares = layer.getFeatures()
    for compare in compares:
        geom1=point.geometry().asPoint()
        geom2=compare.geometry().asPoint()
        dist=d.measureLine(geom1,geom2)
        bear=d.bearing(geom1,geom2)
        metrics.update({dist:bear})
    break
#sort dictionary and extract third value, which is the bearing of vessel in radians
i=1
rad=[]
for j in sorted (metrics):
#    print ((j, metrics[j]), end =" ")
    if i == 3:
        rad.append(metrics[j])
    i+=1

features = layer.getFeatures()
for feature in features:
    # retrieve every feature with its geometry and attributes
    # fetch geometry
    geom = feature.geometry()
    geomSingleType = QgsWkbTypes.isSingleType(geom.wkbType())
    if geom.type() == QgsWkbTypes.PointGeometry:
        ptx=geom.asPoint().x()
        pty=geom.asPoint().y()
        pol = QgsGeometry.fromPolygonXY([[QgsPointXY(ptx+6.625,pty+3.225),QgsPointXY(ptx+5.825,pty+4.025),QgsPointXY(ptx-5.825,pty+4.025),QgsPointXY(ptx-6.625,pty+3.225),QgsPointXY(ptx-6.625,pty-3.225),QgsPointXY(ptx-5.825,pty-4.025),QgsPointXY(ptx+5.825,pty-4.025),QgsPointXY(ptx+6.625,pty-3.225)]])
        centroid = pol.centroid().asPoint()
        pol.rotate(rad[0]*57.295779513-90, centroid)
        #add geometry to layer
        vl = QgsVectorLayer("Polygon?crs=epsg:6348", "temp", "memory")
        pr=vl.dataProvider()
        f = QgsFeature()
        f.setGeometry(pol)
        pr.addFeature(f)
        vl.updateExtents()
        QgsProject.instance().addMapLayer(vl)


texts:
I am using 'queryTask.execute' to get features from a layer. And have change the geometry-paths of the features(polyline), trying to save the changes using 'applyEdits'. This is throwing error 'esri.geometry.normalizeCentralMeridian: 'geometryService' argument is missing'. Updating any attribute is working while path edit is not working.

queryTask.execute(query, function(data){
    if(data.features && data.features.length){
            for(let i in data.features){
                connectedBoxIds = data.features[i].attributes.connectedBoxId.split(',');
                pathIndex = Number(connectedBoxIds.indexOf(componentId)) + 1;
                data.features[i].geometry.paths[0][pathIndex] = newVertex;
            }
            var resultLayer = wireLayer.applyEdits(null, data.features, null);
            resultLayer.then(_ => {
                console.log(_);
            }).otherwise(err => console.warn('error', err));

        }
    });
    queryTask.on('error', lang.hitch(this, 'queryTaskErrorHandler'));


texts:
I have a set of municipality maps, with an ID for each municipality, and a Multipolygon for the geometry. In PostGIS, I've put these data into tables. Table 1 has the municipalities in year 2019 (each row is 1 muni with ID and geometry), and Table 2 has the municipalities in year 2005. Spatial indices for both tables. About 400 rows in each table.
Here are some stats on this data:
All in the same projection: WGS 84

NPoints | Avg   | Max    | Min
2019:   | 14378 | 379875 | 9
2005:   | 24675 | 556334 | 5

Num Geometries | Avg   | Max    | Min
2019:          |  63   | 4402   | 1
2005:          | 131   | 9480   | 1

ST_Valid: 4 invalid geometries out of about 400 in each set.
Basically what I want to do is look at the intersection and see how similar the two geographies are across time. So, I wrote the query below.
SELECT m1.munid, ST_Intersection(m1.wkb_geometry, m2.wkb_geometry)
INTO intersect_table
FROM m2019 m1
LEFT JOIN m2005 m2 on m1.munid = m2.munid

Essentially I join based on the muni IDs (which don't change across time), then create the intersection. But this query basically doesn't finish, I've left it for 30 minutes and I can't imagine it should take longer. Am I being unreasonable or am I going about this incorrectly?


texts: How sensitive is Landsat to atmospheric correction?

texts: How important is it to correct Landsat for atmosphere, haze and convert DN to reflectance? What impact does it have on the results, say for many of the vegetation indices?


texts: Calculating several cells into one cell in same table using ArcPy

texts:
I'm very new to Python and scripting in ArcGIS.
I need to calculate some cells from one field in table into another one cell in another field in same table.

Here's what i need:

Here's what i try and it doesnt work.

      A = arcpy.da.SearchCursor("S2_merge2", ["Name1"] , 'FID = 0')
         B = arcpy.da.SearchCursor("S2_merge2", ["Name1"] , 'FID = 1')
         C = arcpy.da.SearchCursor("S2_merge2", ["Name1"] , 'FID = 2')

      cursor = arcpy.da.UpdateCursor("S2_merge2", ["Name5"] , FID = 0')

for row in cursor:
    row[0] = A+B+C
    cursor.updateRow(row)


texts: Optimisation of MySQL JOIN based on a Within match

texts:
I have a table of points and boundaries, and am trying to add the boundary ID which a point is within. However, using a LEFT JOIN ON Within(location, boundary) it is taking about 3.5 hours match 450,000 points against 350 boundaries. Is there a way to optimise this join?

In more detail:

I have two tables in MySQL 5.6, one of which contains points, each stored as a point and the other containing boundaries, each stored as a geometry:

-- Table of locations, around 0.5 million points
CREATE TABLE locations (
  id INT(11) NOT NULL PRIMARY KEY,
  longitude float(11,6) DEFAULT NULL,
  latitude float(10,6) DEFAULT NULL,
  lonLat point NOT NULL DEFAULT '',
  boundaryId INT(11) DEFAULT NULL
) ENGINE=MyISAM;

-- Populate the lonLat field
UPDATE locations SET lonLat = POINTFROMTEXT(CONCAT('point(',longitude,' ',latitude,')')) WHERE longitude IS NOT NULL AND latitude IS NOT NULL;

-- Add spatial index on lonLat
ALTER TABLE locations ADD SPATIAL INDEX(lonLat);

-- Table of around 350 exact boundaries, some overlapping
CREATE TABLE IF NOT EXISTS boundaries (
  id INT(11) NOT NULL PRIMARY KEY,
  llgeom geometry NOT NULL
);

-- Add spatial index on boundary llgeom:
ALTER TABLE `boundaries` ADD SPATIAL(`llgeom`);

I have a query which updates the location table with the boundary ID that the point for that row is within:

UPDATE locations
LEFT JOIN boundaries ON Within(lonLat, llgeom)
SET boundaryId = boundaries.id;

Note that both lonLat and llgeom both have spatial indexes on them already.

With around 450,000 points and 350 geometries, running on MySQL 5.6, this takes around 3.5 hours. Doing a test limited to just 14 rows takes about 2.1 seconds.

If I run an EXPLAIN, this shows that no indexing is being used:

mysql> explain UPDATE locations LEFT JOIN boundaries ON Within(lonLat, llgeom) SET boundaryId = boundaries.id;
| id | select_type | table      | type | possible_keys | key  | key_len | ref  | rows   | Extra                                              |
|  1 | SIMPLE      | locations  | ALL  | NULL          | NULL | NULL    | NULL | 451010 | NULL                                               |
|  1 | SIMPLE      | boundaries | ALL  | NULL          | NULL | NULL    | NULL |    353 | Using where; Using join buffer (Block Nested Loop) |

This shows that type is ALL in both cases, which is "the worst join type and usually indicates the lack of appropriate indexes on the table."

Is there some improvement I can make which will give much better performance, using indexes?

NB Using the ST_Within function (which gives true boundaries, rather than simplified bounding-box matching) for those same 14 rows takes much longer, 83 seconds:

UPDATE locations
LEFT JOIN boundaries ON ST_Within(lonLat, llgeom)
SET boundaryId = boundaries.id;

However, I have a routine called reallyWithin which has the same result but takes about 2.3 seconds. But whichever of the three functions is used (the bbox Within, the procedure reallyWithin, or the official ST_Within), this works out too slow for 450,000 points.


texts: Splitting landcover polygons (parcels) into equal sizes, while retaining road access for all?

texts:
I would like to split landcover polygons which surround a road polygon into equal sizes, while ALL resulting polygons should have access to the road!

Is there any chance to do this automatically? All programs, licenses etc. are welcome. I just don't have ANY idea how to solve this problem.

To make it more clear, I have attached a basic sketch.


texts: Deleting the "cached" field names for a large CSV file

texts:
Using ArcGIS Desktop Advanced 10.1

I have a large CSV file (36MB with about 600,000 records) that contains single addresses that I want to geocode.  The file is called "SingleAddresses.csv"  The file contains 5 fields.  When I first loaded the file into ArcMap, I had misnamed some of the field names, and so my field names were "SegID", "FromStNo", "Zone1", "Zone2", "Zone3".   So I closed ArcMap, reopened my file in Notepad++, changed the field names to "Address", "City", "Zone1", "Zone2",  "Zone3".  However, when I open up ArcMap again and add the SingleAddresses.csv file again, the old (original) field names show up.

It seems as though ArcMap is caching the field names for the "C:\My Files\SingleAddresses.csv" file.  If I move that file to a different directory and then open it in ArcMap, the field names show up correctly.  If I keep that file in the same directory but rename it (SingleAddresses2.csv) and then open it in ArcMap, the field names show up correctly.  But as long as I keep the same file name in the same directory, it still shows the old field names, even if I delete the original file and create a new empty SingleAddresses.csv file.  It also still shows the old field names if I import that file into a new (different) map.

Note that I cannot duplicate this if the original file is a small file.  The problem only occurs if the file that I originally use is "big"...

Is there a cache somewhere that I can delete that would eliminate this problem?  It seems as though ArcMap caches the field names for a hard-coded file name.  It is really annoying...


texts: Check if button is pressed in QGIS 3 module

texts:
There is a plugin I wrote for QGIS and I added some buttons there.
From main launch Python script I need to check is specific button is clicked and according to that do some stuff.

I saw a default condition in main script:

def run(self):
    if self.first_start == True:
        self.first_start = False
        self.dlg = frame_builderDialog()
    self.dlg.show()
    result = self.dlg.exec_()
    if result:
        pass
        'or do something'

In this case I can catch event if I had an "OK\Cancel" button. But I have another simple buttons. I tried to catch them using something like

result_btn = self.dlg.pushButton.clicked

where pushButton is one of the buttons taken from script with dialog, but nothing were returned.

What should be done here?


texts:
I am new to using GIS and struggling to insert one of my layers that I uploaded in GeoServer in an OpenLayers 3 map. Unfortunately, I do not get the layer displayed.
Can you help me fix the problem?
See attached my code:
var map = new ol.Map({
  target: 'map',
  layers: [
    new ol.layer.Image ({
      source: new ol.source.ImageWMS ({
        url: 'http://localhost:8080/geoserver/wms',
        params: {
          'LAYERS':'cite:SWE_adm2'
        }
      }),
      ratio: 2,
      opacity: 1,
      serverType: 'geoserver'
    })
  ],
  view: new ol.View({
    center: ([55.11, 13.82]),
    zoom: 4,
  })
});


texts: How to request GeoServer elevation data from GeoServer to PNG file with the minimum loss of precision

texts:
I'm requesting alt data on GeoServer from a layer group as a PNG file. The colormap of the style associated with the layer is the following :
<ColorMap>
     <ColorMapEntry color="#000000" quantity="-27" label="values" />// altitude min
     <ColorMapEntry color="#FF0000" quantity="4787"/> // altitude max
</ColorMap>

Which means the altitude is interpolated from black for min alt -27m to red for max alt 4787m.
As it is encoded in a 8bit png file, I only have 256 value between red and black and as consequences for my PNG, I only have a vertical precision of (4787 - (-27))/256 = 18 meters which is not quite good for what I want.
Are there other technique to have a better vertical precision in my PNG file to tackle this problem and take full advantage of the 256 x 256 x 256 = 16 777 216 values of the RGB channel which could theorically give me a much better vertical precision.


texts: Unable to load .gpx file from Python Console

texts:
I'm trying to create a plugin which will extract the gpx files and load them onto the active layer window in QGIS. While I'm able to load it using the QGIS GUI from Layer>New Vector Layer, I'm not able to make it happen from the python console.

You can find a screenshot of the error here.

I've tried using the following commands but in vain.

>>> qgis.utils.iface.addVectorLayer(path,'layer_name','gpx')

But,

>>> vlayer = QgsVectorLayer(path[0],'new_layer','gpx')

did not display any error at the same time it did not appear in the list of active layers. Tried using setActiveLayer() but it returned False.

Any help or suggestions would be welcomed !

Thanks


texts: Visualizing signal strength of mobile towers

texts:
I have a list of latitude and longitude along with signal strength obtained from mobile towers in a certain region. I have all the data in CSV format. In total, I have 4 CSV files. I have mapped the data from CSV file.
I have attached the type of result I wish to have. Can anybody let me know if its possible to get such result in QGIS?


texts: Create Featurecollection from GeoJSON containing list of features


texts: POSTGIS workflow for cleaning up a noisy dataset

texts:
I have a conservation dataset that im trying to simplify for the purpose of using it as a GML layer in openlayers, ie reducing the total shape and node count to quite small.

The data has issues: eg-

about 15K orig polygons of which i only need the 500 largest post dissolve.
overnoded, way too detailed
lots of long thin slithers either inside or outside the main polygons which im not that interested in
hundreds of invalid geometrys, rings with less than three points, self intersections etc

Using this as an exercise whats a good way to clean and simplify a dataset like this?

Ive been trying variations on this so far:

st_simplifypreservetopology,10, to reduce time taken below
hand fix invalids
buffering +20m to swallow the internal slithers,
st_union, trys ot make single feature, mem exhausts usually
unbuffering 40m, to swallow external slithers.
st_union,
buffer +20 to return to right area.
st_simplify,10
st_dump
bin any with area less than 1000ha

Its slow and simplify always creates yet more invalid shapes.

A subquestion is how to make postgis do the equivalent of qgis 'buffer and dissolve', ie merge shapes that adjoin or overlap, but leave other shapes as seperate features.


texts: Which map projection should I use for a country land cover map?

texts:
I am currently completing a piece of GIS work and I am a bit stuck, even after doing some reading.

What would the appropriate map projection and co-ordinate system be for a land cover map for a country the size of Kenya, for example?

Also, what would be the appropriate scale for the digital map data sets that would be used in this application?


texts: Install a specific nightly version of QGIS 3 on Ubuntu?

texts:
Is it possible to install the nightly version of QGIS 3 for a specific date (code version)?

I'm working with QGIS 3 development version which recently has some issues with projections and crashes. I was able to recover a virtual machine (Ubuntu 16.04) in Nov that still works.

Before redoing everything else on the VM, I'd like to ask:

Is it possible to tell apt-get to install a specific nightly version (e.g. for Nov. 21)?

Note: I can't go back to stable versions (QGIS 2) now as I have written some new python 3 code that only works under qgis3.

Related:
installing older version of QGIS on Ubuntu


texts: Create Auto-Incremented Field within SDE Feature Class (SQL Server)

texts:
Has anyone implemented a way of creating an auto-incremented integer field within an ArcSDE (SQL Server 2008 R2) feature class / table?

I don't want to rely on the OBJECTID field because if a feature class is exported to another geodatabase, the OBJECTIDs get re-calculated.  I'm looking for a permanently stored integer field that increments as records are added to the feature class / table.


texts:
I am trying to import a shapefile to my database in postgres but there is no shapefile loader appearing in the plugin menu in PgAdmin.

Postgis is definitely installed as I have added it as an extension to my database

The postgis extension has already been created.  But I continue to get a syntax error:

postgres=# CREATE EXTENSION postgis;
ERROR:  extension "postgis" already exists
postgres=#
postgres=# shp2pgsql -I -s 4326 /Users/alexander/Downloads/pluto/pluto.shp pluto | psql -U postgres -d postgres;
ERROR:  syntax error at or near "shp2pgsql"
LINE 1: shp2pgsql -I -s 4326 /Users/alexander/Downloads/pluto/pluto....
        ^
postgres=#

Any idea as to what the problem could be?


texts: Displaying large number of polygons on web map?

texts:
I've started working on a project with a huge number of polygons, and I'm unsure of how to display/host it. The map I'd like to end up with will have three layers:

1) Shaded polygons representing each of the 70K+ U.S. census tracts. I currently have this as both a shapefile and a spqllite file in Tilemill.

2) Another 15K points representing locations on roadways (this layer will be interactive)

3) A street layer, such as OpenStreetMap

I would set the opacity of the first layer so that the street layer remains visible. I've looked at possibly using UTFGrids, but haven't found much documentation.

Any ideas on how I might want to implement this one?

UPDATE: I'd prefer to use TileMill, but am open to different platforms.


texts: Errors while importing QGIS 3.10 from standalone Python, but not from QGIS's Python Console

texts:
I try to load a QGIS project in standalone Python on Windows 10. This is the code that I use:
from qgis.core import *

qgs = QgsApplication([], False)
QgsApplication.setPrefixPath("C:\Program Files\QGIS 3.10\apps\qgis", True)
QgsApplication.initQgis()
project = QgsProject.instance()
project.read("filename.qgz")
layers = project.mapLayers()
print(layers)

Running this code from the Python console inside the QGIS application works perfectly and it returns a dict of the layers in the project.
However my goal is to build a standalone Python script. I use the Python interpreter that came with the QGIS installation (C:\Program Files\QGIS 3.10\bin\python3.exe) for this. The environment variables are set correctly using e.g. a modification of this description on how to set them. I also double checked that the variables inside the QGIS application and in the command line are the same. Now when executing the above code with the standalone Python interpreter, I get this error message:
Logged warning: Loading a file that was saved with an older version of qgis (saved in 3.4.8-Madeira, loaded in 3.10.2-A Coru?a). Problems may occur.
pj_open_lib(proj.db): call fopen(C:/Users/ESPPC25/AppData/Roaming/python3\profiles\default/proj\proj.db) - succeeded
libpng warning: iCCP: known incorrect sRGB profile
libpng warning: iCCP: known incorrect sRGB profile
QObject::connect: invalid null parameter
qt.svg: Cannot open file ':/images/composer/missing_image.svg', because: No such file or directory
qt.svg: Cannot open file ':/images/composer/missing_image.svg', because: No such file or directory
qt.svg: Cannot open file ':/images/composer/missing_image.svg', because: No such file or directory
qt.svg: Cannot open file ':/images/composer/missing_image.svg', because: No such file or directory
qt.svg: Cannot open file ':/images/composer/missing_image.svg', because: No such file or directory
qt.svg: Cannot open file ':/images/composer/missing_image.svg', because: No such file or directory
qt.svg: Cannot open file ':/images/composer/missing_image.svg', because: No such file or directory
qt.svg: Cannot open file ':/images/composer/missing_image.svg', because: No such file or directory
qt.svg: Cannot open file ':/images/composer/missing_image.svg', because: No such file or directory
QThreadStorage: Thread 0x21e690e9c30 exited after QThreadStorage 4 destroyed
QWaitCondition: Destroyed while threads are still waiting

I can still execute methods like layers.keys(), put printing anything from the layer dict results in the above error message.
What is going on? How can I get the layers from my project in standalone Python? What else could be different between the Python console inside the QGIS application and the standalone Python interpreter?


texts: ArcGIS composite locator results out of order

texts:
I have a composite address locator that cascades through three sub-locators. The issue is when I run the Geocode Addresses tool, it chooses the best match based on the order of locators, not the score of the candidates. So let's say I have these results:

Locator 1: score of 90%
Locator 2: score of 100%
Locator 3: score of 75%

It will always choose the match from Locator 1 (provided it's above the minimum score threshold), even though there's a better fit from Locator 2. Similarly, if I publish this as an ArcGIS Server REST service the results are ordered by locator. This seems like a bug to me -- does anyone know of a workaround? Is this configurable somewhere?


texts: makefile generated using QGIS Plugin builder cannot find pyuic4

texts:
I have QGIS 2.4 installed using the OSGeo4W installer. I have the bin folder in my path and I can run pyuic4 successfully on the command line using osgeo4w.bat, but when I try to run the makefile generated by the QGIS plugin builder, I get an error saying:

pyuic4 -o ui_LoginDialog.py ui_LoginDialog.ui
make: pyuic4: Command not found
make: *** [ui_LoginDialog.py] Error 127

I know the syntax is correct because I can run the first line of that command manually in the command line (under osgeo4w.bat).

Is it possible that the items I have in my PATH are conflicting with the osgeo4w environment? I don't understand how I can run it manually and then it fails in the makefile, in the same osgeo4w session.

My makefile:

UIFILES := $(wildcard *.ui)
PYFILES := $(UIFILES:.ui=.py)

.PHONY: all
all: $(PYFILES)

%.py: %.ui
    pyuic4 -o $@ $<


texts: How do I stop event propagation with rightclick on Leaflet marker?

texts:
I didn't find a rightclick event in Leaflet markers so I started listening to the mousedown event. In the method that resolves this event I check

if (event.originalEvent.button == 2) {
    do_something();
    // now i want to stop propagation of the rightclick
    return;
}

This does work in that do_something() is called, but the event propagates to the map and a contextmenu is being shown. How do I prevent the event to propagate to the map?

Does this work in all major (i.e. IE) browsers?


texts: How to find which areas/regions a list of lat/lons are in?

texts:
The Australian Bureau of Statistics (ABS) issues Esri shapefiles or MapInfo files which have a code # for each geographical area at various 'resolutions' and I have a database with a column of variables and their lat/lon (it's a .csv file). How is it best to get a list of the #'d geographical areas (i.e. the ABS ID_codes) that each variable lat/lon point is within? And how can I get these both in the same database?

I am using QGIS.


texts: Does Google Earth Builder support server-side geoprocessing?

texts:
Is there a feature in Google Earth Builder that supports server-side geoprocessing?

For example, say I upload two large geo datasets to Google's cloud: one with lots of points and another with lots of polylines.

Is there some way I can write a script and upload it to the cloud that:

finds all points in dataset1 within
some distance of a line in dataset2
having some attribute = some
specified value.
find all points within some distance
of the points from step 1 where some
attribute = some value
buffer the points from 2 based on
some other attribute, and return the
results

I would like for this script to run on the server - not on the client.

Update: I watched most of the webinar at Directions Magazine, last week.  The  first 3/4's of the webinar, which included the entire demo and didn't address this question.   Right now it is in the process of being archived, but should be available for viewing again shortly.

Update  The webinar archive is now available for viewing.


texts: Reading CRS and Affine transformation from a NetCDF with CF conventions

texts:
I have a NetCDF with geographical information. I'd like to obtain a GeoTIFF file with rasterio. The problem is how to translate the geographical information from CF standard to something that GeoTIFF can understand.
In the netCDF I have this variable:
int projection_lambert ;
    projection_lambert:grid_mapping_name = "lambert_conformal_conic" ;
    projection_lambert:latitude_of_projection_origin = 63. ;
    projection_lambert:longitude_of_central_meridian = 15. ;
    projection_lambert:standard_parallel = 63., 63. ;
    projection_lambert:proj4 = "+proj=lcc +lat_0=63 +lon_0=15 +lat_1=63 +lat_2=63 +no_defs +R=6.371e+06" ;

Therefore I can create the CRS with:
crs = CRS.from_proj4("+proj=lcc +lat_0=63 +lon_0=15 +lat_1=63 +lat_2=63 +no_defs +R=6.371e+06")

OK, I know the projection, but what is missing is the Affine transformation. This information is implicitly stored in the lat and lon variables the NetCDF also includes, but I'm not sure how to go from there, to the Affine transformation.


texts:
I have two shapefiles, a polygon shapefile (lake) and a line shapefile (boudary)
I want to create a new shapefile "polygon to the north of the boundary"

What tool can I use in QGIS or other free GIS for this task?


texts: How to plot a shapefile over a raster in Python

texts:
I am beginner to work with GIS, so be patient.
My first task is to do:

Read a raster file   (Done)
Create a shape file  (work on it)
Create a new raster from the shapefile (not started)

I am having trouble with to plot the shapefile that I have created over the raster image. The bounding box do not appear at image. How can I figure this out?
I was trying to do something similar to this post: Plot shapefile on top of raster using plot and imshow from matplotlib
from collections import OrderedDict
import rasterio as rio # funcao similar ao R do raster
import rasterio.plot as show
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.collections as cplt
import fiona
import shapely
from descartes import PolygonPatch
from osgeo import gdal

# Read the raster file
filename = "C:/Users/adolf/Desktop/Annia_teste/ALF_16JUN26141551_M2AS_R02C1_pansharpen.tif"
dataset = rio.open(filename)

# Get the band, width and height from raster data
nband   = dataset.count
nwidth  = dataset.width
nheight = dataset.height

# O tipo de dado fica na propriedade .dtypes

# =============================================================================
# Dataset georeferencing #

# Item 1 -  Obtendo o valor corresponde aos vertices da imagem em metros

bounds = dataset.bounds    # Get the spatial geopoint
matrix = dataset.transform # Get the affine matriz to relate the pixel with spatial position
crs    = dataset.crs       # Get Coordinate reference system

# =============================================================================
# Using Fiona to create the shapefile #

# Define your schema as a polygon geom with a couple of fields
schema = {
    'geometry': 'Polygon',
    'properties': OrderedDict([
        ('somecol', 'str'),
        ('anothercol', 'int')
  ])
}

# Create a set of shapefile which are construct by defining a 256 x 256 pixel

square_size = 256
nrow = int(nheight/square_size)
ncol = int(nwidth/square_size)

nrow = 1 #bypass test
ncol = 1

with fiona.open('oo.shp', 'w', driver='ESRI Shapefile',  crs=crs, schema=schema) as c:

    id = 0
    for col_id in range(0,ncol):
        for row_id in range(0,nrow):

            #Criando o poligo a partir dos vertices
            right_bottom = matrix*(row_id*square_size,col_id*square_size)
            left_bottom  = matrix*(row_id*square_size,(col_id+1)*square_size)
            right_top    = matrix*((row_id+1)*square_size,col_id*square_size)
            left_top     = matrix*((row_id+1)*square_size,(col_id+1)*square_size)

            # Corrigindo a ordem p/ escrever os pontos na sequencia correta
            #lr, ur, ll, ul = zip(*polygon)
            polygon = [(right_bottom, right_top, left_top, left_bottom)]

            record = {
                'geometry': {'coordinates': polygon, 'type': 'Polygon'},
                'id': id,
                'properties': OrderedDict([('somecol', 'Some Value'),
                                           ('anothercol', 12345)
                                           ]),
                'type': 'Feature'}
            #Escrevendo no artigo
            c.write(record)
            id = id + 1

# Reading the shape file I have just created
with fiona.open("oo.shp", "r") as shapefile:
    features = [feature["geometry"] for feature in shapefile]

patches = [PolygonPatch(feature) for feature in features]

ax = rio.plot.show(dataset, adjust='linear')
ax.add_collection(cplt.PatchCollection(patches))


texts: Why isn't my ArcGIS JavaScript QueryTask working?

texts:
var map;
    require([
            "esri/map",
            "esri/layers/ArcGISDynamicMapServiceLayer",
            "esri/layers/ImageParameters",
            "esri/tasks/QueryTask",
            "esri/tasks/query",
            "esri/symbols/SimpleMarkerSymbol",
            "esri/InfoTemplate",
            "dojo/_base/Color",
            "dojo/domReady!"
            ], function(
            Map,
            ArcGISDynamicMapServiceLayer,
            ImageParameters,
            QueryTask,
            Query,
            SimpleMarkerSymbol,
            InfoTemplate,
            Color
            ) {

    map = new Map("mapDiv", {
        center: [-123.8425, 45.4552],
        zoom: 16,
        basemap: "streets"
        });

    var imageParameters = new ImageParameters();

    imageParameters.layerOption = ImageParameters.LAYER_OPTION_SHOW;
    imageParameters.layerIds = [4];
    imageParameters.transparent = true;

    var dynamicMapServiceLayer = new ArcGISDynamicMapServiceLayer("http://host:6080/arcgis/rest/services/Data/MapServer", {
          "opacity" : 1.0,
          "imageParameters": imageParameters
        });

    map.addLayer(dynamicMapServiceLayer);

    queryTask = new QueryTask("http://host:6080/arcgis/rest/services/Data/MapServer/4");

    query = new Query();
    query.returnGeometry = true;
    query.outFields = ["OBJECTID"];

    infoTemplate = new InfoTemplate("${OBJECTID}", "OBJECTID : ${OBJECTID}");

    symbol = new SimpleMarkerSymbol();
    symbol.setStyle(SimpleMarkerSymbol.STYLE_SQUARE);
    symbol.setSize(1000);
    symbol.setColor(new Color([255,0,0,1]));

    query.where = "OBJECTID = 1";

    queryTask.execute(query, showResults);

    function showResults(featureSet) {
        //remove all graphics on the maps graphics layer
        map.graphics.clear();

        //Performance enhancer - assign featureSet array to a single variable.
        var resultFeatures = featureSet.features;

        //Loop through each feature returned
        for (var i=0, il=resultFeatures.length; i<il; i++) {
            //Get the current feature from the featureSet.
            //Feature is a graphic
            var graphic = resultFeatures[i];
            graphic.setSymbol(symbol);

            //Set the infoTemplate.
            graphic.setInfoTemplate(infoTemplate);

            //Add graphic to the map graphics layer.
            map.graphics.add(graphic);
        }
    }

});

I followed this ESRI example:

https://developers.arcgis.com/javascript/jshelp/intro_querytask.html

But when I test it I get an error:

Resource interpreted as Script but transferred with MIME type text/plain: "http://host:6080/arcgis/rest/services/Data/MapServer/4/q…Fields=OBJECTID&callback=dojo.io.script.jsonp_dojoIoScript3._jsonpCallback". init.js:495
Resource interpreted as Script but transferred with MIME type text/plain: "http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer?f=json&callback=dojo.io.script.jsonp_dojoIoScript1._jsonpCallback". init.js:495
TypeError {stack: (...), message: "Cannot read property 'clear' of null"}
 "TypeError: Cannot read property 'clear' of null
    at showResults (http://test/test/script.js:61:15)
    at e._successHandler (http://js.arcgis.com/3.10/init.js:490:193)
    at e._handler (http://js.arcgis.com/3.10/init.js:1424:385)
    at http://js.arcgis.com/3.10/init.js:174:23
    at Object.d.load (http://js.arcgis.com/3.10/init.js:1420:411)
    at http://js.arcgis.com/3.10/init.js:630:478
    at c (http://js.arcgis.com/3.10/init.js:74:221)
    at d (http://js.arcgis.com/3.10/init.js:74:10)
    at resolve.callback (http://js.arcgis.com/3.10/init.js:75:350)
    at c (http://js.arcgis.com/3.10/init.js:74:436)"

I know there should be a query result because I tested it through REST. I entered the query OBJECTID = 1 through REST and it returned a single record.

Any idea what could be wrong? I'm stumped.


texts: Open QGIS project on different computers

texts:
I have multiple QGIS projects that need to go to different computers. I put all the files on a USB and loaded it under "Documents" on the other computers like it is saved on my own computer.

However, when opening the projects on another computer the "Handle bad layers" window pops up. There are hundreds of files for each project.

I am wondering if there is a way to get around this without having to click on each file individually and finding the new location?


texts: R- alphahull - duplicate points error help

texts:
I came across R function -alphahull- and believe it is the best option to turn gps points into a polygon (3 or 4 days worth of gps points across residence/neighborhood for each of 300 subjects).
However, I am just starting to use R and I would appreciate any help with removing duplicate points, as I am getting the following error: "Error in tri.mesh(X) : duplicate data points" after running:

library(sp)
library(alphahull)
library(maptools)

#define projection
project2<-"+proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# read point shapefile
data.shape <- readShapePoints("D:/clayera/UWO/msthesistake2/tesisGIS/D1001_GPS_WalkExtSel.shp",proj4string=CRS(project2))
plot(data.shape)
summary(data.shape)
# extract coordinates
x.coords <- coordinates(data.shape) #this way we get duplicate coordinates
# alpha-shape: 100 meter threshold
x.as <- ashape(x.coords[,1], x.coords[,2], alpha=100)


texts: Google Earth Engine - Get Image with highest max Scene NDVI from a time period

texts:
I have a shapefile with several fields, corresponding to different wheat fields of a certain year. Using as Input an image collection converted to NDVIs, I would like to calculate the mean NDVI value for each field for every available date. Then, export the raster of the maximum "mean value" of each field.
The key here is not to do a pixel-wise maximum NDVI selection, but by FIELD.

In other words, I need to get the Sentinel-2 NDVI image from a field that corresponds to the highest mean NDVI value (of the field polygon) registered within a time period. This is not the max() function over a ImageCollection (which synthesizes a image with pixels from different dates). I need the image of a REAL scene, which actually happened at one DATE.

var geometry = /* color: #d63000 */ee.FeatureCollection(
    [ee.Feature(
        ee.Geometry.Polygon(
            [[[-61.58584513408002, -35.77369642524422],
              [-61.58035197001752, -35.778222714936696],
              [-61.5727130387431, -35.772373307284425],
              [-61.56704821330365, -35.77683003785912],
              [-61.55262865763959, -35.76506094262748],
              [-61.5643016312724, -35.75586730183619],
              [-61.58472933512982, -35.77286077435719]]]),
        {
          "system:index": "0"
        }),
    ee.Feature(
        ee.Geometry.Polygon(
            [[[-60.64656168176598, -35.72030669380027],
              [-60.62716394617028, -35.70358079489495],
              [-60.645703374881215, -35.68963986472292],
              [-60.652913152713246, -35.69214941201044],
              [-60.66733270837731, -35.70427777742714]]]),
        {
          "system:index": "1"
        })]);

Map.addLayer(geometry,{},'ROIs')

function getNDVI10 (image){
  return image.normalizedDifference(['B8','B4'])
              .select([0],['NDVI'])
              .addBands(image.metadata('system:time_start').rename('DATE'))
}
function getStats (image){
  return image.select(['NDVI']).reduceRegions({collection: geometry,reducer: ee.Reducer.mean(),scale: 10,
})}

var nbrParams = {min: 0.4, max: 0.8, 'palette': ['red','yellow','green']};
var s2 = ee.ImageCollection("COPERNICUS/S2")
  .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 5))
  .filterBounds(geometry);
var NDVI10_2017 = s2.filterDate('2017-10-01','2017-12-01').map(getNDVI10)

////////////////////////////////////////////////////////////////////////////////////////////////////////////
var NDVIpixel_2017 = NDVI10_2017.qualityMosaic('NDVI').clip(geometry)
Map.addLayer(NDVIpixel_2017.select(['NDVI']),{}, "maxNDVI S2",true)
////////////////////////////////////////////////////////////////////////////////////////////////////////////

var NDVI10_2017_stats = NDVI10_2017.map(getStats)


texts: Caching / Preloading data for use with WMS queries

texts:
Hoping someone can help me with this issue. I've got a map implemented in my web application using OpenLayers that makes WMS queries to GeoServer that in turn queries a Postgres (PostGIS) database to produce WMS tiles. It also makes use of WMS GetFeatureInfo to get information about features being moused over.

Because of the amount of data that is being queried in the initial load for the tiles can take anywhere between 1-10+ minutes depending on the data that is being queried (query can be adjusted by date/time ranges and a few other properties), which is acceptable because it's a lot of data. But I've got a couple of performance related issues that I'm trying to get around.

Every time the zoom level is changed the database has to be re-queried. Is there a way the data can be cached so that it doesn’t have to be re-queried for each change in zoom level? The data doesn’t change depending on zoom level so the query being run is returning the exact same result each time.
When a user mouses over a feature, a GetFeatureInfo request is made and a popup is shown over the map detailing that feature. It once again has to run the whole query again and then it narrows it down to the feature that has been moused over.

Previously before changing to using WMS we generated a single KML file from the data we got from our database query and then this would be used for loading the map at all zoom levels and getting the feature info. We had to change from this approach to deal with bigger and more complex polygons.

So simply, is there anyway to get GeoServer to use a cache result of the query data or is there any way I could perhaps generate the KML file how I was before but then pass it to GeoServer to render it as WMS tiles?

Would much appreciate any input on the matter or any suggestions to get around the problems I've stated.


texts: Arcpy - (Error 999999) Table Already Exists Project Management

texts:
My code is supposed to clip features within a temporary directory, project them and place them into a specific feature dataset within a geodatabase. Then the code goes back and deletes the temporary files in the temporary directory.

The code runs through a list to do this. The first iteration through works flawlessly but when the second iteration begins I receive the following error:

Traceback (most recent call last):
  File "Q:\Admin\AOI_Generator\Scripts\Create_AOI_v1.1_Pick.py", line 219, in
    arcpy.Project_management(temp + '\' + fc + "_" + ds + ".shp", data + "\" + ds + "\" + fc, outCS)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\management.py", line 7822, in Project
    raise e
ExecuteError: ERROR 999999: Error executing function.
The table already exists.
Failed to execute (Project).

I already know that the error occurs during the project_management, I am just unsure why it happens.

Please see the code below.

for buff, ds in n_list:
    def listFcsInGDB(gdb):
            env.workspace = gdb
            fcs = []
            for fds in arcpy.ListDatasets('','feature') + ['']:
                    for fc in arcpy.ListFeatureClasses('','',fds):
                            fcs.append(os.path.join(fds,fc))
            return fcs
    fcs = listFcsInGDB(gdb)
    for fc in fcs:

            arcpy.Clip_analysis(fc,r"Q:\\Admin\\AOI_Generator\\Shapefiles\\Buffer\\" + buff, temp + '\\' + fc + "_" + ds)
            arcpy.Project_management(temp + '\\' + fc + "_" + ds + ".shp", data + "\\" + ds + "\\" + fc, outCS)
            arcpy.Delete_management(temp + '\\' + fc + "_" + ds + ".shp")
    del fcs


texts: Comparing subtype and domain functionality?

texts:
I was trying to determine the best way to standardize my attribute table.

I thought I would use a coded value domain but I am unable to do that without creating a subtype first.  I don't understand this but figured I would just create a blank subtype to create the coded value domain. The trouble is that subtypes require a Long integer field which, when exported to a shapefile, are meaningless, just a number in a field.

Consequently, I decided to just create a new subtype properly. I attempted to delete the before mentioned domain and I get the error:

  Failed to delete a domain from the database. The domain is used by an
  attribute rule.

I removed that "attribute rule" but it still will not delete.

The idea of subtypes and domains are a great idea but this locking problem and the subtype field requirement of being an integer rather than String is an issue for me.

For now I will just create the subtype and leave the domain.


texts: Getting Spreadsheet Data of Points from PNG using ArcMap?

texts:
I'm rather new to using ArcMap.

For my research project, I have a layer of point data. I have been downloading datasets and layers of aggregated annual precipitation across the United States.

All I need to accomplish is to simply get a spreadsheet file (.xls, .csv, or similar formats) that simply contains the rainfall (or the range of rainfall) that point falls into on the map.

So, in the end, all I would need is a spreadsheet that has data that looks like this:
point0 | 0-100
point1 | 251-500
point2 | 1,501-2,000
point3 | 0-100
...

Can anyone help me with this?

The precipitation layer in the image above is a .PNG file.


texts: Using Calculated Fields in ArcGIS Pro Reports

texts:
I have built out a report in ArcGIS Pro 2.4.1 that displays fields from one attribute table. This was a report that I had built out in ArcMap, but we are working on moving everything over to Pro. The issue I am having is that in ArcMap you were able to perform calculations directly on the report utilizing fields.
Ex: =field1/field2
In Pro I have not been able to find a way to do this, so I have had to write a python script to add the field to the table, calculate the field, run the report, and finally delete the field from the table. This is painfully slow compared to running the report in ArcMap. I was hoping someone would have an idea of how I can perform a calculation directly on a Pro Report without having to actually add the calculation to my attribute table.


texts:
I want to clip a raster (Lidar LC81140632013197LGN00 which I downloaded from Earth Explorer) with  shapefile in QGIS version 2.0.1.

Which plugin in QGIS could I use it?


texts: Joining data points to matching polygons

texts: I have polygon data on coral reefs at great barrier reef and separate point data for areas of coral bleaching. I have created a clip that only includes points that exist with in the reef polygons. However I need to be able to identify polygons that have bleaching points in them from those that do not. So for example I would like to be able to click on a polygon that has a bleaching point in it and it automatically select all the other polygons that do?


texts: Points along one or more lines at a time, with rotation

texts:
I'm using ArcMap 10.3, and I have several lines, that I would like to add arrows along at a defined interval, with that lines rotation, so the arrows point correct way along the line.

When I use the construct point (while having a line marked) I can get points at a defined interval for one line at a time, but it does not rotate according to the line direction. The funny thing is, it pretty much shows a preview, when you use this tool, with exactly the type of points I want (maybe a different arrow type), question is, how do I get it as shown on the picture below?

So to sum up, I would like to add points (in the shape of arrows) along a line, with the arrows pointing in the line direction.

Tried what is suggested at:
Unable to rotate Marker Line Symbols

However this gives the difference, that smaller lines will have perhaps 5 markers on a line of 100 meters, but a line of 10000 meters will also have 5, so the interval is not ~ close to each other.


texts: Deleting input shapefile after being used in QGIS modeler

texts:
I have a shapefile that is being used as an input into the start of a model. The model is very simple but makes up a much larger model.

The input is reprojected in the first step but is not needed ever again. Is there a way I can simply delete that shapefile permanently from my computer within the modeler itself?


texts: Is there any packages to download landsat image by shapefile boundary in R or python?

texts: I want to find a package which can automatically download images within a predefined administration boundary (shapefile format). Does anyone know a package can conduct or simplify the process?


texts: Confusing result of raster calculator in QGIS

texts:
I'm confused about the results that I get from the raster calculator in QGIS.

I'm trying to delete the cells of a raster who are prone to flooding. However, a part of my raster is wrongly classified as a high risk area, thus for this part I want to use the original raster and for the other areas I want to delete the cells with a high flood risk.

Below, I added my rasterized polygone layer, red is the part that I want to keep from the original layer and has a value of -2, the rest is -1. This layer is my original layer:

and this is the raster adapted for flood risk:

White areas are areas without data.

In the calculator I typed:

("rasterized" = -2) * "original layer" + ("rasterized" = -1) * "adapted layer"

However, it returns the adapted layer only. Also when I switch the values, I get the adapted layer.

When I make two seperate layers: one with the part of the original layer that I want and one with the part of the adapted layer, this formulas seems to work, but again when I add them I only get the adapted layer.

Does anyone has an idea of what I might be doing wrong?


texts: Peculiar behavior of CRS Selector dialog in QGIS 3.0

texts:
I installed QGIS 3.0 on my computer and noted that when I accessed the CRS Selector dialog by any method it was very slow, and it seemed to slow my whole system (in any case, during the time I was waiting for the CRS dialog to appear, my mouse wouldn't move or it moved every few seconds in large steps).

Two installations on other computers did not have this issue, and the CRS Selector dialog opened immediately. I completely removed QGIS 3.0 and QGIS 2.18 (using Revo Unistaller), rebooted, and reinstalled QGIS 3.0. I got the same behavior, so it seems to me something on my system not related to QGIS is the culprit but I am at a loss as to what it may be.


texts: Automating selection and sliver elimination using QGIS

texts:
At the moment I do this manually:

open a shapefile
open the attribute table
use "select features using an expression" - the expression is "area" < 200

then I use vector/geoprocessing tools/eliminate sliver polygons
I select "common boundary", add a name and location for a new shapefile.

I have tried to automate this process using the modeller, but even after looking at some tutorials I could not achieve this.

Since I have many shapefiles, I would like to process them one at the time automatically.
Note that the resulting shapefile should have the original shapefile name with a prefix or a suffix - I don't mind which.

After spending some time with QGIS, I had a look at the "Processing Toolbox" and I have tried "Eliminate Sliver Polygon" as a Batch Process.
It seems like it should do what I am looking for - but I guess that my parameters are not correct and reading the help file did not help me.
There are three entries that I have to fill in:

Selection Attribute
Comparison
Value

In "selection attribute" I have tried "AREA" (with and without quotes), but is this the correct parameter?


texts: PyGRASS r.mapcalc crashes with ERROR: Unable to make mapset element

texts:
I'm trying to load a GeoTIFF and then do some calculations on it with r.mapcalc

the code is as follows:

from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules import Module

geotiff = 'C:\\Users\\ac\\Downloads\\gtfis\\11.tif'
output = 'C:\\Users\\ac\\Documents\\GIS_DataBase\\test.tif'
curr_raster='layer_1'
r.external(input=geotiff, output=curr_raster, overwrite=True)
g.region(raster=curr_raster)
r.mapcalc('layer_1 = if(layer_1 > 0, layer_1, null())', overwrite = True)

r.out_gdal(input = curr_raster, output = output, format = 'GTiff', overwrite = True)

now to me this looks valid but it chrases at the mapcalc step with the following error:
ERROR: Unable to make mapset element C:\Users\ac\Documents\GIS_DataBase\test.tif (C:\Users\ac\Documents\grassdata\newLocation\ac/C:\Users\ac\Documents\GIS_DataBase\test.tif): Invalid argument

and I don't understand why I'm using win10 and grass 7.8


texts: After I use an algorithm like "sink removal" my raster is shifted around 50 km to the south and west! Before this operation I needed to reproject the DEM from MGI / Austria GK West to WGS 84. Does anybody know the problem and can provide a solution?


texts: How to create a saved, consistently available value/color key when styling a choropleth map in QGIS?

texts:
I've created a simple choropleth map representing birthplace data for the US, by PUMAs, for a particular shapefile (e.g., the layer shapefile_1):

I've categorized the value for each of the PUMAs using QGIS's automatic classifier:

Nothing crazy so far, right? But now, I'd like to apply the exact same palette that QGIS created for shapefile_1, and use it for shapefile_2.

Is this possible?

Can I somehow export this value/color pairing and use it in multiple projects?

Do I need to have both shapefile_1 and shapefile_2 loaded in simultaneously, and somehow style/classify both at the same time?

Is there another solution?

I've struggled with this on multiple projects.


texts: How to displaying u and v components of wind data from a grib file as geoserver coverage?

texts:
I have generated a coverage view from a grib file (u,v components), but I can't style the result of both components.

Is there any way to style both wind components (u and v)?

My style is the next one:

<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>Many color gradient</Name>
<UserStyle>
<Title>SLD Cook Book: Many color gradient</Title>
<FeatureTypeStyle>
<Rule>
<RasterSymbolizer>
   <ChannelSelection>
            <GrayChannel>
                    <SourceChannelName>2</SourceChannelName>
            </GrayChannel>
    </ChannelSelection>
<ColorMap type="ramp" extended="true">
<ColorMapEntry color="#6470a9" quantity="0"/>
<ColorMapEntry color="#5979a9" quantity="2"/>
<ColorMapEntry color="#5997a6" quantity="6"/>
<ColorMapEntry color="#629683" quantity="10"/>
<ColorMapEntry color="#59a159" quantity="14"/>
<ColorMapEntry color="#5ba759" quantity="17"/>
<ColorMapEntry color="#a89e59" quantity="21"/>
<ColorMapEntry color="#a88e59" quantity="25"/>
<ColorMapEntry color="#9d6f61" quantity="29"/>
<ColorMapEntry color="#915968" quantity="33"/>
<ColorMapEntry color="#a75986" quantity="37"/>
<ColorMapEntry color="#85689e" quantity="41"/>
<ColorMapEntry color="#72699e" quantity="47"/>
<ColorMapEntry color="#667f9f" quantity="52"/>
<ColorMapEntry color="#6b949e" quantity="56"/>
<ColorMapEntry color="#8659a9" quantity="194"/>
</ColorMap>
</RasterSymbolizer>
</Rule>

</FeatureTypeStyle>
</UserStyle>
</NamedLayer>
</StyledLayerDescriptor>

If i change  the 1 then i can see the other component but i can't mix them.
Any ideas will be very helpful.


texts:
The tiles of the basemap cannot be downloaded due to a security issue
  while the map is not secured.

Despite the fact that the basemap was shared with everybody on the ArcGIS Online (AGOL), the ArcGIS collector is unable to download the map to be used offline and displays the following message “the map failed to download with the message: Unauthorized access to a secure service“. It seems that there is a security issue that prevents the tiles to be downloaded to the device. What might be the source of it.


texts: QGIS. Polygon vertices automatic bind to WFS point features

texts:
I have a layer with polygons in QGIS.
I have a layer in the same project with points from a WFS.

Is it possible to bind vertices from a polygon to points from the WFS layer. So if the admin of the WFS updates a position of a point, the polygon automaticly changes the bound vertices to the new position of the points?