texts: Calculate field using a field from another feature class

texts:
I have two feature classes, both of them are line features. I'm trying to the traffic count field (AADT) in the Road FC into the Bridges FC based on the distance between the two FCs, this distance being 15m. The other condition in assigning the AADT field to the Bridges FC is that road with the largest AADT field that is within the 15m threshold is chosen. I've already came up with code that assigns a score to each feature in the Bridges FC based on these AADT values (while still following the conditions I mentioned) below:

#Tuples that are used to calculate the score based on AADT, a and b are AADT ranges, c is
the score that gets assigned to the bridge
AADT_Tuple = [(0, 3999, 1),
            (4000, 7999, 2),
            (8000, 13999, 3),
            (14000, 19999, 4),
            (20000, 10000000, 5)]

##Tuple Calculation which does the same thing as the tuple above but for AADT counts
for aadts in AADT_Tuple:
    (minaadt, maxaadt, score) = aadts
    arcpy.MakeFeatureLayer_management(roads, "justselected", "\"RM_AADT\" >="+str(minaadt)+ " AND \"RM_AADT\" <="+str(maxaadt))
    ##AADT = long("RM_AADT")
    arcpy.SelectLayerByLocation_management(bridges,"WITHIN_A_DISTANCE","justselected",15,"NEW_SELECTION")
    arcpy.CalculateField_management(bridges, 'Soc_Traffic_Criteria', score)
    ##arcpy.CalculateField_management(bridges, 'Soc_Traffic_Criteria', AADT)
    arcpy.SelectLayerByAttribute_management(bridges, "CLEAR_SELECTION")
    arcpy.Delete_management("justselected")

The two commented out lines in the for loop (the AADT variable and the 2nd calculate field) were my attempts to do what I'm trying to do but it didn't work. I was thinking of using either Update/Search cursors but I'm not sure how to implement it while still keeping the above code's functionality.


texts:
I have a requirement of GIS to be implemented using Java. I have read the developer guide for OpenStreetMap but there is insufficient information for its set-up technically.

I am completely new to GIS. Any help regarding technical set-up of OpenStreetMap as well as information about using them will be very much appreciated.


texts: Are the CORINE Landcover datasets available by country?

texts:
I found the overall dataset for Europe but it is a huge file for my computer to run. I would like to have a shapefile for the Netherlands specifically.

http://land.copernicus.eu/pan-european/corine-land-cover/clc-2012/view

Hope you can help


texts: Generalizing / simplifying multipolygon in PostGIS without missing polygons?

texts:
By following the rules from this site, I got the output but some polygons do not come within the simplified polygon. I think the problem is with topology.

Is there any solution?


texts: How to Create Feature Layers and Search by Location for all Feature Classes in Geodatabase?

texts:
I'm very new to python, but I'm having trouble finding the answer on my own.

I need to loop through a given file geodatabase, select any features that intersect a project polygon and then search the selection for specific information.

Broadly, how does one create feature layers for all the feature classes in a feature dataset and then pass those layers onto select by location?

My code so far ...

import arcpy
import csv
arcpy.env.overwriteOutput = True

arcpy.env.workspace = "C:\\Users\\Katherine.idea-PC\\Desktop\\FinalProject\\Water.gdb\\WaterDistribution"
featureClassList = arcpy.ListFeatureClasses()
projectArea = "C:\\Users\\Katherine.idea-PC\\Desktop\\FinalProject\\ProjectArea.shp"
outputCSV = "C:\\Users\Katherine.idea-PC\\Desktop\\Final Project\\Output.csv"

for featureClass in featureClassList:
     arcpy.MakeFeatureLayer_management(fc, fc+"_lyr")
     arcpy.SelectLayerByLocation_management (fc + "_lyr", "WITHIN", projectArea)


texts: Seeking Fréchet distance implementation (Python) to cluster trajectories by similarity

texts:
I am working on a trajectory analysis project using Python and its data science related libraries.

I am planning to implement Fréchet distance algorithm for trajectory analysis, each trajectory has sequence of (x,y) coordinates along with timestamp, speed, dist between consecutive points, etc.

I am seeking to implement Fréchet distance to compare each trajectory with another and group the similar curves/trajectories together and form different clusters.

I would like to have some suggestions on this.

Also, I am open to other similarity measure methods like "Dynamic Time Warping" or "Longest common subsequence".


texts: OGR2OGR how to update a new column by the length of polyline

texts:
Hi I can update new column with area of a polygon object, but don't know, how to update new column in polyline layer with length of polyline.

This works very well for AREA:

for %%f in (*.tab) do (
   ogrinfo %%f -sql "ALTER TABLE %%~nf DROP COLUMN PLO1"
   ogr2ogr %%~nf_new%%~xf %%f -sql "SELECT *, OGR_GEOM_AREA AS PLO1 FROM %%~nf" -f "MapInfo file"
   del %%~nf.* /Q
   ren %%~nf_new.* %%~nf.*
)

How to update Object LENGTH the same way, please?

I can do it in MapInfo:

Update tblnm Set Plo1 = ObjectLen(obj, "m")

But update with OGR would be more usefull for me in this case.


texts: QGIS project crashes when project CRS is changed and all layers are activated

texts:
My QGIS project has four vector layers, three of them are from the NaturalEarth dataset (countries, regions and ocean layer). The other two are shapefiles for plant locations and plant names. Whenever I enable on-the-fly transformation from wgs 84 to the respective UTM 33/34/35N for the Scandinavian countries, and activate the country and ocean layer, the project crashes and a minidump is created.

The program doesn't crash when only the region layer is activated. I usually set a style filter on the Region layer according to the country I am working on. I am having this problem only for countries like Finland, Sweden and Norway. For other countries, I was able to produce a raster image from the print composer with all layers activated without any issue. I was also able to change the project CRS  to the respective UTM.

I am using the 2.18 standalone version of QGIS.


texts: Why do the Heat Map (QGIS) and Kernel density estimation (SAGA) tools provide different results?

texts:
The Heat Map tool (QGIS) and the Kernel Density Estimation tool (SAGA) provide different results for the same data as shown below.
Why is this happening?
Is one of them more conceptually correct?

Heatmap - QGIS
Radius: 100 m
 - Weight from field: count
 - Output raster size: 1 m
 - Kernel shape: quartic
 - Output value scaling: scaled
 - Fit: I have understood how to nodes

Kernel - SAGA
Radius: 100 m
 - Weight from field: count
 - Output raster size: 1 m
 - Fit: nodes
 - Kernel shape: quartic kernel
 - Out put value scaling: I have understood how to scaled

The difference between the Heatmap tool results (QGIS) and the results of the Kernel D.E. (SAGA) tool can be seen below.


texts: Using makeline() function for more than 2 points in spatialite?

texts:
Is it possible to use makeline() function to make linestring geometries in spatialite for connecting more than 2 sequence points?
Suppose I have four sequence points P1, P2, P3, P4. Or are there any other way to do this instead of using makeline()?


texts: Is it mandatory to specify --style argument and --tag-transform-script in osm2pgsql query for importing .pbf file?

texts:
I am trying to build a tile server that serves tiles of multiple styles(osm and osm-bright). I have followed this link for implementing this. But here I found that osm2pgsql query for importing .pdf file is specifing the --style argument and --tag-transform-script for importing data to postgres database.

osm2pgsql -s -C 300 -c -G --hstore --style openstreetmap-carto.style --tag-
transform-script openstreetmap-carto.lua -d gis -H $HOSTNAME -U postgres
[.osm or .pbf file]

So my question is that - Is it mandatory to use --style argument and --tag-transform-script osm2pgsql? If yes , what should be the alternate way for serving multiple styles?


texts: Reverse direction for arrow keys in the map window in QGIS

texts: Currently the arrow up moves the window up and the data downwards - I would like to reverse this behavior to match the mining software that I am using.  There must be some Python Code somewhere to do this but have not been able to find it.


texts: How can I create GeoTIFFs from OS Terrain 50 data?

texts:
I'm playing around with R Rayshade, which creates lovely 3D terrain images using GeiTIFFs. I've been testing it out with SRTM 1 Arc-second images from USGS, but these are quite low resolution. OS provides large scale elevation data in Terrain 50. These come as ESRI shapefiles; when I load them into QGIS they appear as contour lines. How can I convert these into a DTM TIFF? I have tried the TIN interpolation tool but it crashes QGIS hard every time, so I suspect I'm not doing something right.

The data is also available in GML, ASCII Grid & GML, and Geopackage, if any of those are easier to work with.


texts: I would create a service to geocode from existing shapefiles.  I know that ArcEngine does not allow the creation of Windows Services.  Is there an open source option someone could recommend?


texts: Equal area map projections: Lambert azimuthal or Albers conic?

texts: When making small-scale regional maps on the order of 500-5000 km in extent where an equal-area projection is desirable, which would you prefer, the Lambert azimuthal or Albers conic? I realize this could be a matter of the cartographer's personal preference and graphical constraints (e.g. aspect ratio of the final map), but is there a general rule of thumb about when to choose one projection over another?


texts:
I am currently looking into obtaining additional formal education in the GIS field. I have looked into several university programs and contacted representative to learn more about the details and pricing. There seems to be a fairly large discrepancy in pricing and I am finding it hard to determine if the higher price is worth it. I have searched online for an independent ranking of accredited GIS programs but so far have only come up with http://www.justinholman.com/2016/01/21/2016-rankings-top-10-online-masters-programs-in-gis/. I am unsure as to the credibility of this site.

Is there any professional GIS organization that offers a ranking system for colleges offering GIS certificates or degrees?

If that is not available, then lists of performance statistics such as graduation rates and % of graduates employed in a related field would suffice. This would help me determine perceived value.

I am looking at online programs as relocation is not an option for me.


texts: What kind of methods are most commonly used to interpolate a relief?

texts: I have got a topographic plan with height measurements which I converted to point shape file with height value in attributes. I was wondering, what kind of methods are most commonly used to interpolate a relief? My guess is that for volume measurements it would be TIN interpolation (for example, QGIS plugin Interpolation), but for pretty visualizations, one of the spline interpolators (SAGA from Processing toolbar) might be more appropriate.


texts: qgis: can't find eclipse pydev path using remote debugger in linux enviroment

texts:
I'm trying to configure pydev and eclipse with qgis for remote debug. Everything worked fine in eclipse and i can start the debug server.
 The problem is that Remote debugger plugin requires a pydev path and i can't find it anywhere in a linux enviroment.
I'm using Ubuntu 16.04.4 LTS (Xenial Xerus), Qgis 2.18, eclipse photon and python 2.7.


texts:
I have already read this post but, I am still getting problems with my QGIS installation, and always show error with libgeos-3.4.2.so. I think, is because the default GEOS version installed by yum install osgeo is 3.5.0-1.rhel7.1.

How can I install Install QGIS 3.4.5 LTR on CentOS 7 by using RPM or anything else?

pd: If I try yum install qgis is installed by default 2.14.9-1.el7


texts: select by attribute on shapefiles with open-source python modules

texts: I usually use ArcGis 9.3.1, but when running a python script that includes Select_by_attribute (from gp), the script runs very very slowly. Is there a way I can perform a select-by-attribute on a opensource-python-module (such as shapely?)? Will it be faster?


texts: Supervised Classification in gvsig - fails repeatedly

texts:
I am teaching an ad hoc class in remote sensing, where students are bringing their own laptops

I wanted to go over how to perform a supervised classification in the open source software I had chosen - gvsig, but even the simplest of tasks seem to fail.

I loaded in a Landsat scene, subsetted it to about 10% of the original area, created an ROI of a single small region of about 10 pixels - and tried to perform a supervised classification on it.  It sat on creating spectral signatures for about 10 minutes, and then gave me an error along the lines of the array out of bounds or something similar.

Does anyone know what I might be doing wrong, or alternatively have another suggestion for open source software that can do supervised classifications that will work on both PC AND Mac?

QGIS will with the SCP plug-in, but the plugin appears to be for PC only.


texts: QGIS grabs incorrect field values for categorized style classification

texts:
as seen in the first image, the field values QGIS returns when I attempt a 'select by expression' are values found nowhere in the field in question ('FLD_ZONE'). Instead, it returns a selection of values from an entirely different field ('DFIRM_ID').

I notice a similar thing in the properties panel when I ask the style guide to auto-classify a categorized color scheme based on the values of the FLD_ZONE field (second image). It again returns values from the DFIRM_ID field, only this time there are only 3.

Does anyone have a clue as to why these QGIS interfaces are interpreting values of this field incorrectly??

imgur.com/3A9f0.jpg


texts: Getting 'locate features along route' tool to measure consistently from the start/end of a polyline

texts:
I am having a bit of an issue understanding the 'Locate features along routes' tool in ArcGIS, I am very new to Linear Referencing toolbox.
I have ran an intersection between a line (river network) and polygon shapefile, and produced a points shapefile from this. I would like to work out the distance that each of these points are from the start or end of the polyline that it intersects. I have ran the Locate features along routes tool for this (after converting the river network into a route layer), but the MEAS value outputted seems to inconsistently measure the distance from either the top or the bottom of the polyline. From a bit of reading, this seems to be something to do with the 'M' value. All shapefiles have the same projection and the river network flows in a consistent direction. I have ran this tool with all of the defaults ticked, and with 'M direction offsetting' unticked, and I still get inconsistent outputs whether the MEAS is taken relative to the top or bottom of the line.
Could anyone point out where I might be going wrong?


texts: Out of range error when converting from geometry to geography in PostGIS

texts:
I am attempting to convert a geometry column in the postgres database to a geography column. I created a new column like this:

alter table scene add column footprint_new geography(POLYGON, 4326);

and I'm trying to populate the new column like so:

update scene set footprint_new = footprint::geography;

and I get the following error:

ERROR:  Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type

Now, it isn't too surprising that there are some out of range values in this database. I want to skip those out of range geometries during the update and identify them later so they can be fixed. How can I do that?


texts:
I was wondering if there is a way to split a shapefile's features in QGIS?

I have many polygons in a shapefile which are separated geometrically. But in their attribute table they all have came in just one row. So I can not use "split vector layer" command. Is there any automated command to split these features?

In the picture as you can see all these polygons are in one layer. I want each part to be in a separated layer.

i tried using multipart to single part command but it raised this error:

Traceback (most recent call last): File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\qgis\VectorSplit.py", line 114, in processAlgorithm sink, dest = QgsProcessingUtils.createFeatureSink(fName, context, fields, geomType, crs) _core.QgsProcessingException: Could not create layer C:/Users/esmae/AppData/Local/Temp/processing_0e22d2fe65794057b6026c802fc150a0/63c9824b0ad746159e05fdf243fec6f8/OUTPUT\fid_27.0.gpkg: Creation of field fid failed (OGR error: Wrong field type for fid) –

my qgis version is 3.6


texts: Determining if ArcGIS Server SOM or SOC is installed on machine?

texts:
I need to create a visual map of my current ArcGIS Server system including where software and data are installed. The System comprises of MSSQL, ArcSDE 9.3, ArcIMS, ArcGIS Server 10.0, and various servers which host Server Object Containers, and at least one Server Object Manager.

I can´t see any SOC or SOM in the windows services. Can anyone give me clue how to get this information from my server environment?


texts: ArcGIS Server REST API: at which URL can I acquire a token?

texts:
I'm working with the ArcGIS Server REST API in a python script. In order to access my server, I need to generate a token and put it into the https://myserverurl.com/arcgis/rest/login page, however, I'm unable to find the URL where a token is generated. This page says "To acquire this token, go to https://myserverurl.com/arcgis and enter https://myserverurl.com/arcgis/rest' for the 'Webapp URL' parameter." This link redirects me to https://myserverurl.com/arcgis/sharing/rest/generatetoken, but it is an empty page.

Some URLs I've tried include: https://myserverurl.com/arcgis/rest/generatetoken and https://myserverurl.com/arcgis/rest/authentication, among others. ESRI's documentation asks for http://gisserver.domain.com:6080/arcgis/tokens, but this is not working either.

Am I missing something?


texts: Why can't QGIS handle the Kavrayskiy7 projection although it is supported in proj4?

texts:
I'd like to make a world map using the Kravayskiy7, however QGIS seems to be unable to handle the projection, even though it is supported by proj4. 704900 (as is the EPSG for kav7) is not listed in the CRS and strangely when I try to define it manually using

+proj=kav7 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m

QGIS tells me that it's invalid (site for definition source). I always assumed that QGIS supported all projections that are featured by proj4. So why is QGIS behaving this way and how could I work around this?


texts:
I’m currently working on an iPhone application that shows visitors where to walk in my local town.  Using GPS on the http://www.opencyclemap.org/  map I want the user to be able to see where they are on the trail.  Just like the beeping dot when you use google or apple maps.

What I don’t know is how do I get my app to interface with opencyclemap so that the user can see where they are with the GPS.  What would be the easiest or most appropriate way to achieve this.  I don’t mind if the user is directed out of my app to the open cycle map page.  If they could stay within the app this would be good.

I understand this is a high level question and as such expect a fair amount of questions to drill into specifics.


texts: Get value of a raster layer from actual extent

texts:
I would like to create a plugin where I could use the values of a raster layer from the actual extent, to do some tricks.

For example if I have a raster which is a box [x,y,width,height] : [0,0,1000,1000], I zoom on it to see the AOI so my actual extent begin [50,50,100,100] (for example), how can I use it these values, without reloading the entire image ?

I would like to do really fast calcul (no reloading, no saving, just displaying with vector calcul => Instant on small area).

I can get the value of one pixel doing it like that :

layer = iface.mapCanvas().layers()[0]
extent = iface.mapCanvas().extent()

layer.dataProvider().identify(QgsPoint(extent.xMinimum(), extent.yMinimum()), QgsRaster.IdentifyFormatValue).results()[1]

Is that a way to us QgsRectangle on identify and not QgsPoint ? Because I would like to get the value of the extent rectangle and not from a point only, without doing a for loop of course to fill an array.


texts: when given bearing, angle, and distance, how to find change in location in 3d

texts:
I am helping a drilling team find the surface location representing a point in the ground.  I am given the following values:

surface location latitude and longitude
Drilling depth (length of drill pipe in the ground)
direction (azimuth) of the drill
Inclination vertical deviation from straight down.

so, for example...
if the entry point of a drill is at Lat: 29.45753895, and Long: -096.63472555

and I'm given depth: 1000 ft, at azimuth: 1.45 deg, in direction 45.23 deg

How do I figure out where on the surface does those numbers represent where the drill head is?


texts: Append, Insert or load Data into geodatabase table field

texts:
I am attempting to calculate point stationing (footage) along a pipeline at a set interval via python script,  thing is I'm just now learning to code in Python, I come from a VBA programming background.  In regards to this exercise I will be creating an empty file geodatabase, within this database I will be creating an empty feature class table.  This feature class table will contain an empty field/column named, "Calculated Footage".  My desired result is to populate the content of my feature clas table via Python script such as the ones below.

I have been successful in obtaining the desired results within a python list as shown below but when attempting to do so in a file geodatabase I seem to be at a lost. Simply put how do I go about applying the logic below to the "Calculated Footage" field in my feature class table? In other words how do I populate the "Calculated Footage" field with the result of my while loop?  I'm simply trying to append, load or insert data into the field in question

B_Measure = 0; E_Measure = 5280
Curr_Measure = 0; Interval = 1000
My_List=[]

C_Measure = 0

# Body of the Main program
while abs(E_Measure-C_Measure)>Interval:
    if My_List==[]:
        My_List.append(C_Measure)
        print C_Measure
    else:
        C_Measure=C_Measure + Interval
        My_List.append(C_Measure)
        print C_Measure[/CODE]

I've done this exercise in Excel as well, using the following VBA Code.

Sub Calculate_Footage()
'Declaring variables
Dim X, Y, Beginning_Point, Ending_Point, Interval, Current_Point As Integer

'Assigning initial value to variables
 X = 2
 Y = 1
 Beginning_Point = 0
 Ending_Point = 5280
 Interval = 1000
 Current_Point = 0

'Assign Field/Column Name
Cells(1, 1).Value = "Calculated Footage"

'This While Loop will populate my "Calculated Footage" field
Do While (Ending_Point - Current_Point) > Interval
Current_Point = Current_Point + Interval
Cells(X, Y).Value = Current_Point

'Move to the next row but same column
X = X + 1
Loop

End Sub

Current Python Code
    # Import system modules
    import arcpy
    from arcpy import env

# Set workspace in which to create geodatabase
env.workspace = "C:/data"

# Set local variables
out_folder_path = "C:/data"
out_name = "fGDB.gdb"

# Execute CreateFileGDB
arcpy.CreateFileGDB_management(out_folder_path, out_name)

# Set environment settings for Table in geodatabase
in_database = out_folder_path+"/"+out_name
table_Name = "Custom_Stationing"
arcpy.CreateFeatureclass_management(in_database, table_Name)

# Set workspace for Add Field
in_Table = in_database+"/"+table_Name

Field_Names=["Route_Id","Route_Code","Series","Measure","Calculated_Field"]
for n in Field_Names:
    arcpy.AddField_management(in_Table, n)

I would like to apply the following code sample to insert rows located hereINSert Cursor


texts: Optimal way to estimate date range of an ImageCollection in Earth Engine

texts:
What is a most straightfoward way to get a date range of an ee.ImageCollection as formatted date strings?

For example, ee.ImageCollection('LANDSAT/LT04/C01/T1_32DAY_EVI').get('date_range') returns [398044800000,756259200000], so first I need to cast both items into ee.Date and then format() or print() them to get 1982-08-13 and 1993-12-19, but I think there should be easier ways to estimate date range of a collection.


texts: How to get Lon/Lat values from the postgis raster type column

texts:
I have a tiff raster in my PostGIS 2.0 database. I want to get the value lat/lon and extra data like temperature from the raster type column in postgis table.

I want the table with columns lon,lat,temperature etc.. from raster type column(rast).
  For each corresponding lat and lon threre is a temperature value.

The tiff file details by gdalinfo as:

Driver: GTiff/GeoTIFF
Files: atm_ncep_gfs_temp_20130422_003.tif
atm_ncep_gfs_temp_20130422_003.tif.aux.xml
Size is 720, 361
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
         AUTHORITY["EPSG","7030"]],
       AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (-0.250000000000000,90.250000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
AREA_OR_POINT=Area
latitude#long_name=latitude
latitude#units=degrees_north
longitude#long_name=longitude
longitude#units=degrees_east
NC_GLOBAL#Conventions=COARDS
NC_GLOBAL#GRIB2_grid_template=0
NC_GLOBAL#History=created by wgrib2
time#_FillValue=9.999e+020
time#long_name=verification time generated by wgrib2 function verftime()
time#reference_date=2013.04.22 00:00:00 UTC
time#reference_time=1366588800
time#reference_time_description=forecast or accumulated, reference date is fixed
time#reference_time_type=3
time#time_step=0
time#time_step_setting=auto
time#units=seconds since 1970-01-01 00:00:00.0 0:00
TMP_2maboveground#_FillValue=9.9990003e+020
TMP_2maboveground#level=2 m above ground
TMP_2maboveground#long_name=Temperature
TMP_2maboveground#short_name=TMP_2maboveground
TMP_2maboveground#units=K
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -0.2500000,  90.2500000) (  0d15' 0.00"W, 90d15' 0.00"N)
Lower Left  (  -0.2500000, -90.2500000) (  0d15' 0.00"W, 90d15' 0.00"S)
Upper Right (     359.750,      90.250) (359d45' 0.00"E, 90d15' 0.00"N)
Lower Right (     359.750,     -90.250) (359d45' 0.00"E, 90d15' 0.00"S)
Center      ( 179.7500000,   0.0000000) (179d45' 0.00"E,  0d 0' 0.01"N)
Band 1 Block=720x2 Type=Float32, ColorInterp=Gray
 Min=211.200 Max=311.600
 Minimum=211.200, Maximum=311.600, Mean=277.878, StdDev=21.574
 NoData Value=9.9990002605540088e+020
 Metadata:
_FillValue=9.9990003e+020
 level=2 m above ground
 long_name=Temperature
NETCDF_DIMENSION_time=1366599600
NETCDF_time_units=seconds since 1970-01-01 00:00:00.0 0:00
NETCDF_VARNAME=TMP_2maboveground
short_name=TMP_2maboveground
STATISTICS_MAXIMUM=311.60000610352
STATISTICS_MEAN=277.87782601544
STATISTICS_MINIMUM=211.19999694824
STATISTICS_STDDEV=21.573645315608
units=K

I have:

postgresql-9.2
postgis-2.0 raster support
raster imported into PostGIS with SRID 4326

How can i do this? Thanks in advance!


texts: What Open Source Packages are available to provide routing?

texts:
I am looking for ways to validate mileage charged by our contractors.  They are paid from a central location to a customer, and then from the customer to a drop off location.

I use QGIS for a bunch of stuff and also have access to MapInfo Pro.

I need to do the following:

  route individual trips to get total distance in each direction
  also need to know if parts of one trip overlap the other, and by how much (i.e. km)
  and if the second trip ends anywhere along the route of the first trip (i.e. they came back the same way they went out)

I have been messing around with PostgreSQL and pgRouting for a couple days and it's pretty confusing to me.


texts: Script exported from ModelBuilder gives ERROR 000733 at SelectByLocation?

texts:
I created a model in ModeBbuilder, and I exported the model as a script, but when I try to run it I get the following error. I have added my error below and the code below that. Thanks for any help.

  Error: Traceback (most recent call last):   File
  "C:\Temp\VirtualCampus\IGIS\IGISprogrammingscript_1.py", line 43, in

      arcpy.MakeFeatureLayer_management(Points_of_interest, Pointsofinterest_Layer, "", "", "FID FID VISIBLE NONE;Shape Shape
  VISIBLE NONE;CATEGORY CATEGORY VISIBLE NONE;NAME NAME VISIBLE NONE")
  File "c:\program files\arcgis\desktop10.1\arcpy\arcpy\management.py",
  line 5774, in MakeFeatureLayer
      raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000733: Output Layer: Same as input
  C:\Temp\VirtualCampus\IGIS\Super Output Areas.shp Failed to execute
  (MakeFeatureLayer).

Here is my code:

~# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# IGISprogrammingscript_1.py
# Created on: 2013-07-09 21:39:46.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: IGISprogrammingscript_1 <Thames_shp> <Points_of_interest> <Super_Output_Areas_shp> <Thames_Buffer_shp> <superoutputareas_clip_shp>
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Script arguments
Thames_shp = arcpy.GetParameterAsText(0)
if Thames_shp == '#' or not Thames_shp:
    Thames_shp = "C:\\Temp\\VirtualCampus\\IGIS\\Thames.shp" # provide a default value if unspecified

Points_of_interest = arcpy.GetParameterAsText(1)
if Points_of_interest == '#' or not Points_of_interest:
    Points_of_interest = "Points of interest" # provide a default value if unspecified

Super_Output_Areas_shp = arcpy.GetParameterAsText(2)
if Super_Output_Areas_shp == '#' or not Super_Output_Areas_shp:
    Super_Output_Areas_shp = "C:\\Temp\\VirtualCampus\\IGIS\\Super Output Areas.shp" # provide a default value if unspecified

Thames_Buffer_shp = arcpy.GetParameterAsText(3)
if Thames_Buffer_shp == '#' or not Thames_Buffer_shp:
    Thames_Buffer_shp = "C:\\Temp\\VirtualCampus\\IGIS\\Thames_Buffer.shp" # provide a default value if unspecified

superoutputareas_clip_shp = arcpy.GetParameterAsText(4)
if superoutputareas_clip_shp == '#' or not superoutputareas_clip_shp:
    superoutputareas_clip_shp = "C:\\Temp\\VirtualCampus\\IGIS\\superoutputareas_clip.shp" # provide a default value if unspecified

# Local variables:
Pointsofinterest_Layer__2_ = Thames_Buffer_shp
poi_cat = Pointsofinterest_Layer__2_
poi_cat_View = poi_cat
Pointsofinterest_Layer = Points_of_interest
pop = superoutputareas_clip_shp
pop_View = pop

# Process: Make Feature Layer
arcpy.MakeFeatureLayer_management(Points_of_interest, Pointsofinterest_Layer, "", "", "FID FID VISIBLE NONE;Shape Shape VISIBLE NONE;CATEGORY CATEGORY VISIBLE NONE;NAME NAME VISIBLE NONE")

# Process: Buffer
arcpy.Buffer_analysis(Thames_shp, Thames_Buffer_shp, "300 Meters", "FULL", "ROUND", "NONE", "")

# Process: Select Layer By Location
arcpy.SelectLayerByLocation_management(Pointsofinterest_Layer, "COMPLETELY_WITHIN", Thames_Buffer_shp, "", "NEW_SELECTION")

# Process: Summary Statistics
arcpy.Statistics_analysis(Pointsofinterest_Layer__2_, poi_cat, "CATEGORY COUNT", "CATEGORY")

# Process: Make Table View
arcpy.MakeTableView_management(poi_cat, poi_cat_View, "", "", "CATEGORY CATEGORY VISIBLE NONE;FREQUENCY FREQUENCY VISIBLE NONE;COUNT_CATEGORY COUNT_CATEGORY VISIBLE NONE")

# Process: Clip
arcpy.Clip_analysis(Super_Output_Areas_shp, Thames_Buffer_shp, superoutputareas_clip_shp, "")

# Process: Summary Statistics (2)
arcpy.Statistics_analysis(superoutputareas_clip_shp, pop, "LSOA_POP SUM", "")

# Process: Make Table View (2)
arcpy.MakeTableView_management(pop, pop_View, "", "", "Rowid Rowid VISIBLE NONE;FID FID VISIBLE NONE;FREQUENCY FREQUENCY VISIBLE NONE;SUM_LSOA_POP SUM_LSOA_POP VISIBLE NONE")~


texts:
I need to convert a Google Maps polygon object (can has multi paths and holes) to wkt string (multipolygon) at client.  Then i'll able to save polygon as geography in Sql Server 2012.

I tried Wicket (https://github.com/arthur-e/Wicket) but it has a bug converting multipolygons.

Is there any ready to use library?
Any ideas, please?

Thanks.


texts: Using copySelectionToClipboard, pasteFromClipboard in PyQGIS

texts:
I am trying to write a simple script to copy the selected features from the active layer to a new memory layer. This would be equivalent to Edit -> Copy Features, Edit -> Paste Features As -> Temporary Scratch Layer.
The code I have so far is
#Get the type of the active layer
active_layer = iface.activeLayer()
geom_type = active_layer.geometryType()
if geom_type == QgsWkbTypes.PointGeometry:
    type_string = "Point"
elif geom_type == QgsWkbTypes.LineGeometry:
    type_string = "LineString"
elif geom_type == QgsWkbTypes.PolygonGeometry:
    type_string = "Polygon"
else:
    type_string = "unexpected"

#create a url from the layer type and the layer crs
crs = active_layer.sourceCrs();
url = type_string + '?' + crs.authid()

#create the new memory layer
new_layer = QgsVectorLayer(url, "copy", "memory")

#add the attributes of the active layer to the new memory layer`
attrs = active_layer.dataProvider().fields().toList()
new_layer.dataProvider().addAttributes(attrs)
new_layer.updateFields()

#copy the selecton from the active layer and paste to the new memory layer, add to TOC
iface.copySelectionToClipboard(active_layer)
iface.pasteFromClipboard(new_layer)
QgsProject.instance().addMapLayer(new_layer)

When executed the new layer is created, appears in the layer list,  and the attribute table shows the appropriate fields. A pop-up at the top of the QGIS canvas says "n features were successfully pasted", but the attribute table of the new memory layer remains empty.
What key step(s) am I missing?
(This is just an exercise in trying to learn a bit about PyQGIS, so any other comments on approach are more than welcome).


texts:
I downloaded some Landsat 8 OLI/TIRS C1 Level-1 images from EarthExplorer. I was hoping to get reflectance data out of the blue, green and red band; however, when I opened those images in ArcMap, the scale does not seem to be reflectance, given that is not in the range between 0 - 1 (theoretical values of reflectance).

The values that I am getting are 0 - 40079 (blue band) and 0 - 40336 (red band). My question is, is there some type of corrections that I need to do? maybe some transformation?


texts:
I've used QGIS2WEB to output a leaflet map, all of which works well. It's a simple map showing markers, I want the markers to update based on a csv which is updated daily. Omnivore looks like my solution, so i've tried to edit the code QGIS2WEB gives me.

I can see that the CSV in QGIS is converted to a geojson file (i think) and stored in the data folder as;

csv_data.js

This is then called in the main html file using

<script src="data/csv_data.js"></script>.

I've tried editing the contents of the csv_data.js from;

var csv_data = {
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name":
"urn:ogc:def:crs:OGC:1.3:CRS84" } },

to

var csv_data = omnivore.csv(url)

I'm using the mapbox omnivore library with a call to that in the main body of the html using;

<script src='https://api.mapbox.com/mapbox.js/plugins/leaflet-omnivore/v0.3.1/leaflet-omnivore.min.js'></script>

But this returns no markers, the folder structure is as QGIS2WEB output as screen grab.

In the console i get two errors;

Uncaught ReferenceError: omnivore is not defined &
Uncaught ReferenceError: csv_data is not defined

I'm new to javascript and html so not sure whhat is causing either of these


texts: Modifying script published as Geoprocessing Service in ArcGIS for Server?

texts:
In ArcGIS 10.1 before publishing script as a new GPService I have to run it in ArcMap. It is a good solution when I handle small amount of data. When I want to change something in my script  I have to run all processes again. I have found my script in ArcGIS Server in GPServer catalog in C:\arcgisserver\directories\arcgissystem\arcgisinput. When I try to edit it I cannot save any changes.

I tried to solve this problem by stopping service and changing permissions to catalogs. It is annoying to run all processes in order to change one small thing.

Does anyone have idea how to save changes?


texts: Return a set of unique pixel values in a raster with PyQGIS?

texts: I am trying to test whether a value exists in a raster against all the values present.  I am using pyqgis (Python 2.7) in in QGIS 2.14.  I have found a number of ways to get point location values and even to list all values, but I cannot find a way to return a set of unique values.


texts: Raster cells with no data being generated by las2dem?

texts:
Anyone has an idea where these lines (with no elevation data) come from? They are generated while using the las2dem tool from LAStools.

What needs to be done to avoid this? The raw data (.las file) does not contain data with value 0.


texts: Counting nodata and data pixels using PostGIS

texts:
I'm using PostgreSQL 9.2.4 and PostGIS 2.0.3.

I have some raster layer and I'm trying to get the count of all pixels (including nodata), cliped by polygon.

Raster was imported by raster2sql:

raster2pgsql -I -d -C -M -s ESPG_code -t 80x80 input_source > output.sql
psql -h localhost -U user_name -d database -f output.sql

So my code is:

WITH
    polygon As (
        SELECT ST_GeometryFromText('wtk_polygon_geometry))', epsg) As geom
    )
SELECT
    ST_Count(ST_Clip(r.rast,p.geom),1,false) As val
FROM
    raster As r,
    polygon As p
WHERE
    ST_Intersects(r.rast,p.geom)

and I'm geting overestimated value and it looks like below.

but the query should output a count value looks like this - non highlighted pixels:


texts: Redefining nodata values into zero in QGIS without changing the other values?

texts:
I have a raster layer with nodata pixels and I want these to be 0. I have tried the SAGA tool "reclassify values", the GDAL "fill nodata" (also with a validity mask containing only zeroes) and the SAGA "close gaps".

My problem is that with all these, the other values are also being changed --for example, the maximum value in my initial raster layer is 32.94 but after the reclassification, the maximum value is 29.4.
I don't understand why the other values are also being transformed!?


texts:
In postgres, I created a view with a sequence nextval() to have an unique id (called vid).
I tried to vizualise this view in QGIS but it failed. It seems that nexval() causes some issue according to a message :

  "error : can't execute nextval() in a read only transaction"

How this issue can be solved please?


texts: Where to find .mxd examples for ArcGIS 10

texts:
I'm trying to learn how to use ArcMap.

I would like to find some example MXD files for ArcMap, but couldn't find anything googling.

Does anyone know where I can find any MXD samples for ArcGIS 10?


texts: Find closest point to shapefile coastline in Python

texts:
I have a xarray (674 lats & 488 Lons) and want to find the closest distance for each point to the coastline in meters.
I found this solution:
Finding closest point to shapefile coastline Python
which is basically what I want to do. However, the distance is measured in degrees and not in meters (see here).
I could convert degrees to meter by using 1deg=111km but this would not be very accurate for larger domains and domains further south.
My working example is below:
import geopandas as gpd
from shapely.geometry import Point, box
from random import uniform
from concurrent.futures import ThreadPoolExecutor
from tqdm.notebook import tqdm
import cartopy
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import pandas as pd

lon = np.arange(129.4, 153.75+0.05, 0.05)
lat = np.arange(-43.75, -10.1+0.05, 0.05)

precip = 10 * np.random.rand(len(lat), len(lon))

ds = xr.Dataset({"precip": (["lat", "lon"], precip)},coords={"lon": lon,"lat": lat})

ds['precip'].plot()

def get_distance_to_coast(arr):

    def compute_distance(point):
        point['dist_to_coastline'] = point['geometry'].distance(coastline)
        return point

    print('Get shape file...')
    world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

    #single geom for Norway
    aus = world[world["name"]=="Australia"].dissolve(by='name').iloc[0].geometry

    #single geom for the coastline
    c = cartopy.io.shapereader.natural_earth(resolution='50m', category='physical', name='coastline')

    c     = gpd.read_file(c)
    c.crs = 'EPSG:4326'

    print('Get coastline...')
    coastline = gpd.clip(c.to_crs('EPSG:4326'), aus.buffer(0.25)).iloc[0].geometry

    print('Group lat/lon points...')
    points = []
    i = 0
    for ilat in arr['lat']:
        for ilon in arr['lon']:
                points.append({'id':i, 'geometry':Point(ilon,ilat)})
                i+=1

    print('Computing distances...')
    with ThreadPoolExecutor(max_workers=4) as tpe:
        result = list(tqdm(tpe.map(compute_distance, points), desc="computing distances", total=len(points)))

    gdf = gpd.GeoDataFrame.from_records(result)

    print('Convert to xarray...')
    lon = gdf['geometry'].x
    lat = gdf['geometry'].y
    df1 = pd.DataFrame(gdf)
    df1['lat'] = lat
    df1['lon'] = lon
    df1 = df1.drop(columns=['id','geometry'])
    df1 = df1.set_index(['lat', 'lon'])
    xarr = df1.to_xarray()

    return xarr

dist = get_distance_to_coast(ds['precip'])

plt.figure()
dist['dist_to_coastline'].plot()
plt.show()

My guess is to replace the point['geometry'].distance(coastline) with something using the haversine function, but I have no idea how do to this, especially something halfway efficient.


texts:
How can I use the mother of all gps conversion tools, gpsbabel, to create shapefiles? It reads and writes almost everything. Shapefile is conspicuous by it's absence.

A bit of research turns up that shapefile is intentionally not listed as a supported filetype because "Shapefiles are not well defined.  They're more free-form than most of the formats we support and often require additional arguments or modification to the source." -- Robert Lipe, gpsbabel author

However the same thread goes on to say shapefile is "in the '-^3' output and marked 'hidden' for this very reason - it's not practical to treat it as 'just another format' "

Suppose I'm willing to take my chances, how does one get at the hidden formats?
(xcsv is mentioned as well).


texts: Performing Boundary Select and Sum Attributes in QGIS?

texts:
I´m a MapInfo user, I´ve gone over to QGIS and it´s got a lot going for it, except for some "missing" simple features from MapInfo I just can´t do without.

MapInfo Boundary Select - spatial query on the fly to select objects in one layer contained within objects in another layer.
MapInfo Sum Attributes - show sum of numerical attribute fields on the fly for selected objected.

QGIS seems to have no simple solution to these.   I've tried to dig out every related plugin I can find, with no real luck.


texts: I downloaded data from OpenStreetMap with QGIS, and I need to get the gas stations. Can you tell me what labels tell if something is a gas station?


texts: Getting a masked raster with no NAs and sum pixels within overlaying polygon?

texts:
I have done an IDW interpolation in R for precipitation data following the routine of Manny Gimond (https://mgimond.github.io/Spatial/interpolation-in-r.html) to get precipitation inside a catchment.

Here is the routine:

library(gstat) #Use gstat's idw routine
library(sp)    #Used for the spsample function

#Create an empty grid where n is the total number of cells
grd <- as.data.frame(spsample(P, "regular", n=50000))
names(grd) <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd) <- TRUE  # Create SpatialPixel object
fullgrid(grd) <- TRUE  # Create SpatialGrid object

# Add P's projection information to the empty grid
proj4string(grd) <- proj4string(P)

# Interpolate the grid cells using a power value of 2 (idp=2.0)
P.idw <- gstat::idw(Precip_in ~ 1, P, newdata=grd, idp=2.0)

# Convert to raster object then clip to Texas
r <- raster(P.idw)
r.m <- mask(r, W)

The problem is that in the last step (i.e masking the raster with interpolated values (r) to the boundaries of the geometry of interest) a great number of NA values appear that correspond to point lying out of the masking polygon.

What I am looking for is a function to use in conjunction with mask that ensures that in the final masked raster these NA values are removed so that I have raster with no NA value.


texts: Linking Geoportal Map Preview to a Geoserver service?

texts:
I am very new to Geoportal and Geoserver

I have an installation of Geoportal. In the metadata of a resource I can indicate a URL link. I have seen that depending of the linked data the "Preview link" is shown. I also have a Geoserver serving WCS,WMS and WFS services.

What do I need to put in the link URL of a Geoportal resource so when I click preview the request is sent to Geoserver and a shape or grid is rendered in the preview Window? Or if it has to be WMS which image is supported by the preview window?

Please let me know if I am totally confused about linking Geoportal and Geoserver.

Many thanks,
Carlos.


texts: How to create displacement vectors in QGIS using GPS points with different dates?

texts: Is there a way to calculate and visualize displacement vectors in QGIS ? I have DGPS points taken at the same place (ground control points), but their position moved over time as they are located on a landslide.


texts: Code features based on their rank of an attrbute

texts: I would like to have all the records in a table given a rank value based on the values in another field. The caveat is that there are tied rankings. For example in the attached screen cap, EPDO is the data to rank, and RANK_EPDO is field to store the rank. Anyone aware of any tools in ArcGIS that will accomplish this?


texts: Source for shapefiles or geo-referenced data for Madhya Pradesh, India

texts:
I am looking for data sources/shapefiles (demographics, cities, infrastructure, transportation, land cover, geology, etc.) for the state of Madhya Pradesh, India (specifically the Dewas or Khargone districts near Indore) and wondering if anyone has suggestions for good (re: free) sources for data for this region.

I am having difficulty downloading state from the state forest department (http://www.mpforest.org/intranet/gias/) and the MP GIS site doesn't seem to be free (http://www.nchse.org/mpgis/introduction.htm#.ViwBmCsYFSU).

Any suggestions?


texts:
From some other posts I found that i can determine a layers visiblity by something like

#assume lyr_name the layers name to determine visibity for
l = QgsMapLayerRegistry.instance().mapLayersByName(lyr_name)[0]
vis = QgsProject.instance().layerTreeRoot().findLayer(l.id()).isVisible()

Apart from this being put in a method like

def isVis(lyr_name):
    l = QgsMapLayerRegistry.instance().mapLayersByName(lyr_name)[0]
    return QgsProject.instance().layerTreeRoot().findLayer(l.id()).isVisible()

Is this the recommendend way to check if a layer is visible, whichs name is known or is there a more simpler approach?


texts: gdalwarp cutline creates artefacts not quite no data

texts:
I am trying to clip an imagery tile by a shapefile using gdalwarp.

I have tried a number of different commands...

gdalwarp -of GTiff -cutline Coast.shp IN.JPG OUT.tif

gdalwarp -of GTiff -dstnodata -cutline Coast.shp -crop_to_cutline IN.JPG OUT.tif

However both gdalwarp commands result in some pixels that are not nodata.

I do not think I can use -srcnodata  -dstnodata  because the water is different RGB values so cannot specify a correct value.

How can I clip without these errors?


texts: Difference in area calculation between QGIS and R

texts:
I used R to calculate the area of every land polygon of the world (Projection of the source shapefile: WGS84; CRS: EPGS: 4326).

To do that, I loaded the required packages ("sp" and "raster") and I simply used this line of code:

polygon_name$new_field <- area(polygon_name)/1000000

I divided by 1000000 to obtain the area in square kilometers.

However, repeating the same operation with QGIS 3.0.3 for some random polygons, I observed that the result is slightly different. For example, with R the surface of Antarctica is 13,924,270 kmq, whereas using QGIS is 13,761,440 kmq.

What could cause such difference?


texts: Arranging rulebased symbology alphabeticaly

texts:
QGIS version 3.16.4
I'm working with a layer that has a lot of different items in rulebased symbology. But they seem to be ordered randomly. Is it possible to arrange them alphabeticaly?
I'm hoping that with this I can also have them show up alphabeticaly in the legend, now its just a big mess. The legend shows a smaller selection of the items, so I'm also fine with arranging them myself in the legend...but that's also not possible
Here a picture of the legend I want to make alphabeticaly:

And a picture of the symbology page: (not all items are shown)

And here a picture of what the layers look like in regards to each other:


texts: Getting the coordinates of the points inside of a shapefile/Getting the coordinates of the closest point if there are none

texts:
I'm working with an ERA5 (0.25x0.25 degree) temperature dataset for Europe in order to get average temperature in small statistical regions called NUTS 3. So far I've used a shapefile of the NUTS 3 regions to visualize the temperature on an specific day and time.I also plotted a 0.25x0.25 degree grid. The code can be seen below:
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
from mpl_toolkits.basemap import Basemap

ds=xr.load_dataset('Europe_temp_dataset-test.nc')
ds_celsius=ds-273
ds_celsius.t2m.attrs['units']='Celsius'
df=ds_celsius.to_dataframe()
df.reset_index(inplace=True)
df=df.set_index('time')
df.to_csv('Europe_temperature_dataset.csv')

ds_celsius.to_netcdf('Europe_temp_in_celsius_dataset-test4.nc')  #To solve: If file already exists, it gives [Errno 13] error

data=Dataset('C:/Users/sand_jv/Desktop/NetCDF4/Europe_temp_in_celsius_dataset-test4.nc')
lats=data.variables['latitude'][:]
lons=data.variables['longitude'][:]
time=data.variables['time'][:]
temp=data.variables['t2m'][:]

map=Basemap(projection='merc',
            llcrnrlon=-11.0,
            llcrnrlat=28.0,
            urcrnrlon=32.0,
            urcrnrlat=72.0,
            resolution='i'
            )

lon,lat=np.meshgrid(lons,lats)
x,y=map(lon,lat)

c_scheme=map.pcolor(x,y,np.squeeze(temp[0,:,:]),cmap='jet',shading='auto')
    #map.drawcoastlines()
    #map.drawcountries(linewidth=1,color='k')
map.drawmeridians(np.arange(-11,32,0.25))
map.drawparallels(np.arange(28,72,0.25))

map.readshapefile('NUTS_BN_01M_2016_4326','NUTS 3',color='r')

cbar=map.colorbar(c_scheme,location='right', pad='10%')
plt.title('Temperature 2m, 01.01.2019, 00:00')
plt.show()

The results can be seen in the next picture:
As an example, zooming in, the plot looks like this:
As it can be observed in the second picture, there are two cases:

Regions with one or more grid points inside
Regions with no grid points inside (circled in green in the example)

For case 1, I would like to get the coordinates of every grid point contained in each region so I can get an average of the temperature. For case 2, if there are no grid points inside of the region, I would like to find the coordinates of the closest grid point.
Any ideas on how to achieve this?


texts: What algorithm is used by the $geowithin operator of Mongodb

texts:
MongoDB can perform geospatial query such as finding restaurants inside region with the operator $geowithin. For example:

db.places.find(
   {
     loc: {
       $geoWithin: {
          $geometry: {
             type : "Polygon" ,
             coordinates: [ [ [ 0, 0 ], [ 3, 6 ], [ 6, 1 ], [ 0, 0 ] ] ]
          }
       }
     }
   }
)

I would like to know which algorithm is used by the operator $geoWithin of MongoDB. Can someone provide the pseudo-code?


texts: How to invert polygons (or other options for exterior styling)?

texts:
I have a shapefile containing many polygons that I want to use to mask a raster. I have tried gdalwarp -cutline <shapefile> -crop_to_cutline -srcnodata -9999 -dstnodata -9999 <source_raster> <target_raster>, but the results are visually unappealing because of the low resolution of the raster (i.e. jagged edges). It would look better if I could simply place the vector layer over the raster as a masking layer in the map. I would need to see the underlying raster inside the polygons, but not outside. Two options have occurred to me for accomplishing this.

Invert the polygons and save as new shapefile then fill with opaque background color.
Style the vector layer such that inside the polygons is transparent and outside is opaque.

Being new to GIS, I haven't figured out how to accomplish either of these approaches. I have at my disposal the gdal framework, QGIS, and GRASS. I think that I would prefer a GRASS command-line solution, as I may need to repeat this process more than a few thousand times in the next few days. Also, please feel free to direct me to the appropriate documentation that I have obviously missed.

I hope that this question isn't a duplicate of one of the many I've read about masking rasters with vectors. It appears to me that this case is different because I specifically want the final map to display partial pixels from the raster via a smooth cut by a vector.


texts: Style manager not working to batch style group in QGIS 2.0.1

texts:
I am trying to apply batch styling to a group of layers using the style manager by setting the conditions, but the symbology does not seem to be applying. This is considering the multiqml plugin is not ready for QGIS 2.0 yet. In any case, is it possible to also upon adding vector data, specify a preloaded style, i.e. qml style descriptor or sld when importing them into the layer work-space?

For that matter, how do you enable batch labeling for a group?


texts: QGIS 2.6.1 plugin repository, is it broken?

texts:
In QGIS 2.6.1 the plugin repository is unavailable. The URL is:
http://plugins.qgis.org/plugins/plugins.xml
Is there an alternative plugin repository, or a way to resolve it? Have tried two different computers, and still cannot download plugins.


texts: Separate buildings from other objects from a normalized digital surface model

texts: I have created a normalized Digital Surface Model using the substraction DSM-DEM. I have set a threshold value of 2 meters in order to eliminate objects such as shrubs, cars etc. My problem is how to distinguish buildings from the other objects such as trees in an automated way?


texts: Append triangular symbols along polyline in ArcMap?

texts: I have a polyline that needs to have triangles placed along the inside of the line as shown in the sample picture below. I've searched through ArcMap's symbology set and thought I could just add a second line along the original with a symbology of triangle shapes along the new line but cannot find anything appropriate in the symbology set. Any thoughts?


texts:
I've been using the Map My Fitness App from under amour for a bit now and I'd definitely recommend it. Where this app falls short is the export of data. It only goes to a Garmin specific TCX format.

Now I don't know a whole lot about this format outside of its XML-ish. I'm weird and love to collect, store, classify and visualize my runs/hikes/kayaks. I probably spend more time on this then actually exercising :S.

There doesn't seem to be a native ArcGIS import tool or process for getting TCX into a shapefile or feature class. So I've started decoding it and am in the process of writing a parser that with the help of the Shapefile.py module will convert the data.

Now it's taking some time as the TCX XML is well funky. Now I figure I'll probably get the parser done in the next few days when I have some spare time, but has anyone else written any code or built a module or tool within ArcGIS that would accomplish this?


texts: Showing Outline of Island on Topographical Contour Map using QGIS?

texts:
I am fairly new to GIS.

I have been reading through this site and have not quite been able to figure out what I am doing wrong.

I am trying to create a contoured topographical map that I can cut on my laser cutter. The issue is that the area I am working on is coastal and there are lots of areas at or around sea level. As a result, the contour map does not include the general outline of the island/coast and is, consequently, largely unrecognizable.

How do I adjust the settings in QGIS to take into consideration the outline of the land mass I'm working with?


texts: Get XY coordinates in UTM from web-mercator where UTM zone is not known

texts:
How to get XY(northing/easting) coordinates in UTM  where data is WGS_1984_Web_Mercator_Auxiliary_Sphere(WKID: 3857) coordinate system and UTM zone is unknown.
For e.g. India falls under multiple UTM zones and map is in Web_Mercator_Auxiliary_Sphere.
I am trying to replicate Google earth's UTM coordinate display using ArcGIS Javascript API.
Thanks.

Nagesh


texts: When doing a reconcile in ArcGIS ArcMap (v10.2), you get the Reconcile window as shown in the screenshot.  At the bottom it asks how you want to resolve conflicts (in favor of target or edit version).  Once you hit OK on this window, if there are conflicts, you get a conflict resolution window that shows the conflicts and allows you to choose how each conflict can be resolved, so I'm wondering about the purpose of the conflict question in the reconcile window because if there are conflicts, I get the conflict window anyway.  In other words, it seems to be redundant.  I'm sure I'm missing something in the logic here.


texts: Can I change an existing Cartodb map to use a different dataset?

texts: I've created a map with lots of layers and info box info and sql codes, but now I need to use a different dataset from the one I was originally using. Is there a quick way to just swap out the dataset (all the fields are the same)? I haven't been able to find it on any of the documentations on Cartodb.


texts: Export image from Google Earth Engine and then open it in ArcGIS Pro. The pixel size is not correct

texts: I exported an image from Google Earth Engine and then opened it in ArcGIS Pro. The cell size in the layer property was 8.98315284119486E-06. The original cell size should be one meter. The image looks normal in ArcGIS pro. Is this caused by unit difference? How can I fix it?


texts:
I have got the following problem:
I need to change some of the values of one attribute in my point shapefile. I want to do it in modelbuilder, but I don't know how. I have found a script which can do it in Python. The script is as following:

import arcpy
...
... features = arcpy.UpdateCursor(r"C:\Users\peter.valent\Documents\ArcGIS\Data\Vah po LM\test_points.shp")
... for feature in features:
...         feature.precipitation = value
...         features.updateRow(feature)
... del feature,features

This is what I put to modelbuilder

Expression:

 myfcn(r"%file_name%","%value%")

Code Block:

import arcpy

def myfcn(path,value):
features = arcpy.UpdateCursor(path)
for feature in features:
        feature.precipitation = value
        features.updateRow(feature)
del feature,features


texts: gdalwarp: differentiate between empty-in-source and outside-of-cutline?

texts:
I'm using gdalwarp to clip a polygon out of a map layer using --crop_to_cutline.

The original map layer has a lot of empty pixels (well, transparent / alpha=0, to be precise). As it stands, my output is a GeoTIFF, where the pixels that were originally transparent remain as such, and the areas of the image that would fall outside of the cutline are also transparent.

What I want is to be able to differentiate, within that output image, between the pixels that were transparent in the source and the pixels that fell outside of the cutline, say by giving them a specific value that wouldn't exist anywhere else in the image.

Is there an easy way to do that?


texts:
I am wondering if anyone can point me to a resource that will provide a way to model out smells of a garbage dump based on wind? I have been searching through articles on down wind pollution models and the like but have not been able to successfully find anything.

I am fine with using any open source solutions or ESRI based products.


texts:
I have downloaded CSV occurrence datasets for voles and kestrels in the UK, then displayed the XY data of each as points using the latitude and longitude. Not sure if this is the right place to post but I have a few queries

Is there any way to remove certain points? They are very large datasets and I only wish to investigate England and Wales - although I have data for Scotland and Northern Ireland as well. I could only come up with a painstakingly slow method of trawling through 10,000 odd points based on their co-ordinates. Can I somehow delete the XY data points?
Is there a better way to depict the populations other than clumps of dots?
Is there anything cool I can do with the data? (it's for an assignment - we've barely been taught GIS and I would like to get cheeky extra marks) I added locations of SSSI's (highest level of British protected areas) in England and Wales to the maps. Could I do some sort of analysis based upon that?

I am using esri arcgis 10.1


texts: How to convert points distribution toward surface information?

texts:
I'm working with QGIS 2.6

I have a points layer representing different bus and train stops. I applied a categorized symbology; each colour represents a class of value (field regarding time access from a starting point).

I would like represents this time information as different polygons, because in a 2nd time I have to intersect this information with another point layer.

Could you throw a light on this for me? Rasterize this data would be a interesting way? I tried vector -> convex envelope, but the result is not what I am looking for. The spatial distribution of points seems to be too complex to convert the points cover toward convex envelope.

I paste below a SS of my current point data:


texts: How can I save a vector layer to memory (or is there a reason QGIS doesn't allow this)?

texts:
It's impressive how many different vector file formats QGIS (using ver 3.6.0) allows me to save a layer as:

But as far as I can tell, QGIS doesn't allow me to save a vector layer to memory. In my normal workflow in manipulating vector layers, I find myself having to save several intermediate vector layers to file, clogging up my file system or databsase with something like "vector_layer", "vector_layer1", "vector_layer2", "vector_layer2a", etc. It is often trivial to perform the updates/edits to the layer, and because I really don't care about these intermediate files, I'd rather be able to just save them in memory within the QGIS environment. Saving to memory would allow me to simply remove the layers from the project when I finish my workflow instead of having to delete from my Postgres database or my desktop.

When performing a geoprocessing task, there is the option to save to memory, but why can't I have the option to export a layer to memory on its own?


texts: What do the columns in this GIS file mean? Are they distances to 3 reference points?

texts:
I was handed this file and expected to know what to do with it.  Does it possibly mean that point 1 is ~6539 ft from some marker and ~11370 from another marker?  Are these types of file usually in metric or standard?

Sorry for all the sub-questions, I just would like to know basically how to glean info from this file.

Sample Data:

1   6539.603    11370.224   3387.635
2   6558.233    11370.503   3387.023
3   6586.439    11357.517   3386.395
4   6604.131    11357.935   3385.786
5   6615.655    11357.978   3385.310
6   6629.992    11358.072   3385.020
7   6640.550    11358.197   3384.977
8   6659.566    11359.253   3384.952
9   6666.896    11381.068   3385.652
10  6658.798    11381.012   3385.682

LuWi pointed out that the answer would depend on what I need to do with the data:

I am expected to convert this information into latitude and longitude points;  I feel like I need more information from the giver of this data, but I do not know what to ask for.


texts:
I want to get all contiguous polygons for every county. In QGIS I can do this manually with a Spatial Query, but for a large data set this result impractical.

Is there a way to run the spatial query plugin from QGIS python console?


texts: Google MyMaps doesn't show polygons properly

texts:
I have encountered a weird problem. Since I did the map digitizing in QGIS and next exported it to .kml format using MMQGIS plugin I am losing a few shapes in Google MyMaps.
Straight after file export, I opened it on

Google Earth where everything is fine. Unfortunately, when I exported the same file to Google MyMaps, a few polygons are missing. They have been turned into the polylines only.
I want them to be displayed as polygons with proper color fill, set previously yet in QGIS.
I know, that MMQGIS sometimes doesn't return all colors properly, but this time everything seems to be alright with the plugin, as the .kml file is fully valid in Google Earth.
I have 2 questions then:

What should I do to avoid a situation like this? Have I done some mistakes during digitizing?
How can I turn the enclosed polyline into the polygon in Google MyMaps?

Next, when I export the map from Google MyMaps into .kml format and open it in Google Earth, I am getting the same result with lost polygons. I think, that Google MyMaps must have some bug then.


texts: Convert string into date format problems

texts:
I have a string column within an attribute table that has a string in the format of:

yyyyMMdd HH:mm:ss (i.e 20150417 13:32:05)

I've tried using the convert time field tool within the python IDE in ArcMap and it's not giving me the full date/time.

import arcpy as arc

arc.ConvertTimeField_management(in_table="TEST_TIME",
                                input_time_field="Date__Time",
                                input_time_format="yyyyMMdd HH:mm:ss;1033;;",
                                output_time_field="Date_Conv",
                                output_time_type="DATE",
                                output_time_format="yyyyMMdd HH:mm:ss")

I've used both "not used" and the format I want in output_time_format

The resulting output field is always MM/dd/yyyy with no time stamp included.


texts: Would like to install QGIS web app builder

texts:
I would like to install QGIS web app builder plugin. I followed all the instructions in the website as closely as possible (installing Boundless desktop, etc etc.)

When I open QGIS for boundless desktop 1.1 as per this website and try to follow the steps I get the following error when switching to plugins tab in Boundless Connect pane:

  There has been a problem performing the search: Network error #503:
  Error downloading
  https://api.boundlessgeo.com/search/matchAll?version=0.1&cat=BLOG,DOC,QA,LESSON,LC&si=0&c=20
  - server replied:

So instead I tried to install the web app builder using Plugins>Manage and Install Plugins. I did find it in my list of plugins but when I click on Install plugin I get the following error:

  Plugin installation failed:Error downloading
  https://qgis.boundlessgeo.com/plugins/packages-auth/webappbuilder.1.0.zip
  - server replied: FORBIDDEN

There is no material online which indicates why this problem is happening. What is going on here?


texts: Importing raster and vector: GRASS GIS 6.4.4

texts:
Almost giving up on this task, and need some support. I´ve been studying GRASS GIS and training with some tutorials with The NC-dataset. What I want is doing the same task as the one called "Surface and volume calculation" from Open Source GIS - A Grass GIS Approach (3. edition).

My problem is that neither raster or vector formats can be imported. The files are there, but the "import-button" are still greyed out. I´ve got a bigger DTM-file wich i clip within QGIS, plus I need to rasterize my vectorlayer for mask when preparing for r.volume excavation.

What works is the original DTM-file. By that I  mean that it can be imported (consists of an .asc and .tab for mapinfo)

What am I doing wrong?


texts: Get coordinates by mouse click as user's input?

texts:
I'm doing some stuff with digitizing and one thing that would nice to have is to get coordinates by mouse click.

My problem is that I don't know how can I create utility similar to this from map tool, where I can easly create method:

def canvasReleaseEvent(self, e):
    if e.button() == Qt.LeftButton:
        x = e.pos().x()
        y = e.pos().y()

        point = self.iface.mapCanvas().getCoordinateTransform().toMapCoordinates(x, y)
        print(point)

For example, when user wants to digitizing, he has to show a point on canvas, where to start. To do this, he clicks on a button to call a function for getting coordinates from canvas.

The main problem I struggle with is that I don't know how I can craete method for catching coordinates, that only works when I call it (not like in maptool, when it works all the time).

Is there any way to do that?


texts: Creating convex hull with geometry method of ArcPy?

texts:
I am trying to use the geometry method to calculate the convex hull around two polygons (polyA, polyB) with the following code:

# Union both polygons to one (multipart) polygon:
polyAB = polyA.union(polyB)
# Calculating the convex hull around polyAB:
convexhullAB = polyAB.convexHull()

To check if it is working I am writing the combined polygon (polyAB) and the convex hull (convexhullAB) to a shapefile:

cur_poly = arcpy.da.InsertCursor("polyAB.shp", ["SHAPE@"])
cur_poly.insertRow((polyAB,))
cur_convex = arcpy.da.InsertCursor("convex.shp", ["SHAPE@"])
cur_convex.insertRow((convexhullAB,))

The shapefile polyAB.shp contains the (correct) polygon, but convex.shp is empty.
Has anyone an idea what went wrong?


texts: Get rid of jumping GPS data with pykalman

texts:
I have gps data that I get from a smartphone application. Whenever the smartphone is stationary, the gps points are jumping. I understand that the signal is inaccurate due to the reception in a city between buildings and signal loss whenever inside.

From this post I wanted to give a shot to the Kalman filter.  Using this article I was able to try out the Ramer-Douglas-Peucker algorithm on the latitude and longitude, and try the pykalman package for the elevation data. Also I have tried the pykalman example from the github project to play with the filter.

According to these readings, I assume to have the wrong input parameters :

measurements = numpy.column_stack([longitude_list, latitude_list])
# F_k : state transition matrix
F = numpy.array([ [1, 0],
                  [0, 1] ])

# H_k : the observation matrix
H = numpy.array([ [1, 0],
                  [0, 1] ])

# R : covariance R of delta_k observation noise N(0, R)
R = numpy.diag([1e-4, 1e-4])

kf = kf.em(measurements, n_iter=100)
(filtered_state_means, filtered_state_covariances) = kf.filter(measurements)
(smoothed_state_means, smoothed_state_covariances) = kf.smooth(measurements)

The figures below are from matplotlib. The top left is with an iteration km.em(n_iter=2), top right with iteration 10, bottom left is iteration 50, bottom right is iteration 100. Whenever I try higher, I have a timeout. It does not seem that my filter work much on this case. Indeed the same shape is output (it might looks different at first because of the scale axis).

How can I improve my stationary gps jumping data?


texts: what is the stream order coming to next given map

texts:
I am following Strahler's method for stream ordering.

Is the green colored stream 4th order, according to this method?


texts: Sentinel-3 scaling factor in Google Earth Engine

texts:
Does anybody know how to use the scaling factors from Google Earth Engine with Sentinel-3 dataset?
The band values provided in ee.ImageCollection("COPERNICUS/S3/OLCI" by default (GEE, orange) are integers which are very close to the official band values from ESA's Copernicus Open Access Hub (S3, blue) (see figure). The GEE values seem to be the rounded ESA values, thus I assume the scaling factors should bring the decimal part. From my experience the decimal part can also be negative.


texts: Include  in GPX when writing from ogr2ogr

texts:
http://www.gdal.org/drv_gpx.html seems to indicate that ogr2ogr should be able to write <link>[1] elements into a GPX file, however when using ogr2ogr to convert from PostGIS to GPX, I get the message,

Field of name 'link_href' is not supported in GPX schema. Use GPX_USE_EXTENSIONS creation option to allow use of the <extensions> element.

I've tried naming the field, link, href, link_href and link0_href, but none work. Any hints?

[1] http://www.topografix.com/GPX/1/1/#type_linkType

Update: Looking at the sources http://sources.debian.net/src/gdal/1.10.1%2Bdfsg-8/ogr/ogrsf_frmts/gpx/ogrgpxlayer.cpp/ it doesn't seem like it supports writing  elements, only reading, is this correct?


texts: In C# ArcObjects, how to get layer field names from ArcGIS Server layer?

texts:
I'm creating an add-in to ArcMap using C#. I need to get all the fields of every layer.

Code below works to get the fields when the layer is added via a shapefile. However, I can't get the fields when a layer is added from an ArcGIS Server URL (Map service or Feature service), via Catalog -> GIS Servers -> Add ArcGIS Server. The cast of ILayer to ILayerFields comes back null every time. I know the fields exist in the layer, as they show up with an Identify operation.

Is there a different method I should use to get fields from the ArcGIS Server layers? Or is it even possible to get the fields from Server layers?

private void AddLayerNames(string aPath, ILayer aParent, Dictionary<string, List<string>> aFieldNames)
    {
        if (aPath == string.Empty)
        {
            aPath = aParent.Name;
        }
        else
        {
            aPath = aPath + '/' + aParent.Name;
        }

        var lCompositeLayer = aParent as ICompositeLayer;
        if (lCompositeLayer == null)
        {
            if (!aFieldNames.ContainsKey(aPath))
            {
                aFieldNames.Add(aPath, new List<string>());
                var lFields = aParent as ILayerFields;
                if (lFields != null)
                {
                    for (int lFieldIndex = 0; lFieldIndex < lFields.FieldCount; lFieldIndex++)
                    {
                        var lFieldName = lFields.Field[lFieldIndex].Name;
                        aFieldNames[aPath].Add(lFieldName);
                    }
                }

            }
        }
        else
        {
            for (int lLayerIndex = 0; lLayerIndex < lCompositeLayer.Count; lLayerIndex++)
            {
                AddLayerNames(aPath, lCompositeLayer.Layer[lLayerIndex], aFieldNames);
            }
        }
    }


texts: AttributeError: Object: Tool or environment  not found

texts:
I have the following Python toolbox in a file named test.pyt:

import arcpy

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "MyToolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [MyTool]

class MyTool(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "MyTool"
        self.description = "Does some custom things"
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
        return [arcpy.Parameter(
            displayName="My Arg",
            name="my_arg",
            datatype="GPString",
            parameterType="Required",
            direction="Input"
        )]

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

    def updateParameters(self, parameters):
        """Modify the values and properties of parameters before internal
        validation is performed.  This method is called whenever a parameter
        has been changed."""
        return

    def updateMessages(self, parameters):
        """Modify the messages created by internal validation for each tool
        parameter.  This method is called after internal validation."""
        return

    def execute(self, parameters, messages):
        """The source code of the tool."""
        my_arg = parameters[0].valueAsText
        my_arg_str = u'My Arg: {}'.format(my_arg)
        arcpy.AddMessage(my_arg_str)
        print my_arg_str
        return

It works fine when I run it using ArcMap:

But if I try to execute it using a Python script, I get the following error:

Traceback (most recent call last):
  File "C:/Users/MyUser/Documents/ArcGIS/make_package.py", line 13, in <module>
    arcpy.MyTool_mytools(u'some value')
  File "C:\Users\MyUser\Documents\ArcGIS\test.pyt", line 22, in MyTool
    """Define parameter definitions"""
AttributeError: Object: Tool or environment <MyTool_mytools> not found

The following script reproduces the error and shows that arcpy can find the newly imported toolbox:

import inspect
import os

import arcpy

pyt_path = os.path.join(os.path.dirname(__file__), u'test.pyt')
arcpy.ImportToolbox(pyt_path, u'mytools')

print(arcpy.MyTool_mytools)
print(inspect.getargspec(arcpy.MyTool_mytools))

arcpy.MyTool_mytools(u'some value')

print(u'SUCCESS!')

Its output is

<function MyTool at 0x0DBDC7F0>
ArgSpec(args=['my_arg'], varargs=None, keywords=None, defaults=(None,))
Traceback (most recent call last):
  File "C:/Users/Administrator/Documents/ArcGIS/make_package.py", line 13, in <module>
    arcpy.MyTool_mytools(u'some value')
  File "C:\Users\Administrator\Documents\ArcGIS\test.pyt", line 22, in MyTool
    """Define parameter definitions"""
AttributeError: Object: Tool or environment <MyTool_mytools> not found

Note the function and ArgSpec objects. These show that the function exists on the arcpy module. So the failure is somewhere inside those functions.

Why does executing the tool fail?


texts: Missing label using Easy Custom Labeling?

texts:
I’m using the Easy Custom Labeling plugin (with QGIS 1.8). I have one feature label which refused to show itself. It’s just not visible. When I look at the attributes table, the feature name appears in the same columns as the names of all the other features. I’ve tried moving nearby labels away from it--no difference. I’ve tried moving the feature (just to see what would happen)—nothing. If I turn on Labels in the corresponding feature layer, the label appears.

Where’s my label?


texts:
I had run the contour tool in my script to generate polylines.

Contour("rectExtract", "C:/fakepath/Class1.shp", 50,0)

Now I want to label the polylines with the contour-value. If I execute this script:

layer = arcpy.mapping.ListLayers(IMXD, "")[1]
layer.showLabels=True
arcpy.RefreshActiveView()

it will label the polylines with the ID-Value. How can I change this?

Thanks for your help.


texts: GeoServer Image Mosaic Plugin error loading time-series GeoTIFFs

texts:
In an AWS EC2 Instance Linux installed with GeoServer
I am trying to load time-series raster imagery using the Image Mosaic plugin as described here:
https://docs.geoserver.org/latest/en/user/tutorials/imagemosaic_timeseries/imagemosaic_timeseries.html
I have modified the time to -GMT and created all .properties files.
However, when I try to load the directory to publish I am getting the following error

I was checking the read write access of GeoServer to the data folder it  was okay.
Contents of -indexer.properties file:
Caching=false
TimeAttribute=ingestion
ElevationAttribute=elevation
Schema=*the_geom:Polygon,location:String,ingestion:java.util.Date,elevation:Double
PropertyCollectors=TimestampFileNameExtractorSPItimeregex,DoubleFileNameExtractorSPIelevationregex
-timeregx.properties
regex=[0-9]{8}
I do not know where it went wrong?

Yes in my question the indexer.properties  contents are not complete that was a paste error.
I have used the indexer.properties without the elevation collector as suggested. Still, it is not working.
I am attaching the log

Also the filenames I am using in valid form.


texts:
I'm looking at a handful of different mapping projects that I'd like to do something Torque-like with, but I'd like maps to stay on the map for a limited period of time.

Some fake data:

Address                   Start Year         End Year
7706 Heather Limits        2002               2006
5938 Noble Shadow Inlet    2002               2003
9128 Round Pioneer Link    2003               2005
7176 Crystal Via           2004               2005
4475 Harvest Square        2004               2007

Think of a map of buildings participating in a rent subsidy program. One building enters the program in 1988 and exits in 2008. A dozen enter the program in 2002, and will exit in 2032. It seems like I can add points to the map based on a single date, but there's not an easy way to take them off on a specific day, too.

I'd have to turn the "cumulative" option off and then create basically a row for each site for each year that it should be on. Something like

TorqueYear   Address                   Start Year         End Year
2002      7706 Heather Limits        2002               2006
2002      5938 Noble Shadow Inlet    2002               2003
2003      7706 Heather Limits        2002               2006
2003      5938 Noble Shadow Inlet    2002               2003
2003      9128 Round Pioneer Link    2003               2005
2004      7706 Heather Limits        2002               2006
2004      9128 Round Pioneer Link    2003               2005
2004      7176 Crystal Via           2004               2005
2004      4475 Harvest Square        2004               2007
2005      7706 Heather Limits        2002               2006
2005      9128 Round Pioneer Link    2003               2005
2005      7176 Crystal Via           2004               2005
2005      4475 Harvest Square        2004               2007
2006      7706 Heather Limits        2002               2006
2006      4475 Harvest Square        2004               2007
2007 ...

Am I right about that or am I missing something?


texts: Batch process with custom plugins in QGIS

texts:
In QGIS, it is very easy to do batch processing like this but how do you do run a batch process for a custom plugin?

More detail:

I have a number of point shapefiles and want to convert them to lines using the Points2One tool and want to set up a batch process to run the plugin for each of my shapefiles (>100) in a given directory.

Is there a way to do this?


texts: How to select first coordinates in multipoint or in collection used for point reducing?

texts:
I am using PostGis and I have geometry like this:

MULTIPOINT(1673118.56314624 5778263.30559826,1672914.41092247
5778129.131443,1672384.7731665 5778105.44570288,1673513.70724225
5778525.74090022)

OR collection:

GEOMETRYCOLLECTION(POINT(1673118.56314624 5778263.30559826),
POINT(1672914.41092247 5778129.131443),POINT(1672384.7731665 5778105.44570288),
POINT(1673513.70724225 5778525.74090022))

I get it like this for collection:

SELECT st_astext(unnest(ST_ClusterWithin(geom, 700)))
gc FROM table

or like this to get multipoints

SELECT st_astext(ST_CollectionExtract(unnest(ST_ClusterWithin(geom, 700)) , 1))
gc FROM towersnew

And I would like to to transform it to point where the only coordinates are first coordinates (others are omitted). Result is supposed to be:

POINT(1673118.56314624 5778263.30559826)


texts: Have to click twice to open popup leaflet

texts:
I am trying to make a map in which a marker can be clicked and a popup opens up with some information. I have successfully done all of that however my first click on any of the markers result in the popup being placed randomly on the map. If I click again, the popup will place itself correctly (on top of the marker). I believe that I am doing something wrong when calling the "onEachFeature" function. Here is the code:

    // Importing the data from the geojson
$.getJSON("data/factories.geojson", function(data)
{
    factories = L.geoJson(data, {
        pointToLayer: function(feature, latlng)
        {
            return L.circleMarker(latlng, style);
        },
        onEachFeature: onEachFeature,
        style: style
    }).addTo(map);
});
// Action for each feature of the choropleth
function onEachFeature(feature, layer)
{
    console.log("onEachFeature");
    layer.on({
            click: highlightSelection
    });
}
// Highlighting the selected administrative area
function highlightSelection(e)
{
    // Some code her to add text and image to the popup
    factories.bindPopup(image + text);

    **// This will reset the highlight of the previously clicked marker
    if(firstPass == false)
    {
        console.log(lastLayer);
        resetHighlight(feature, lastLayer);
    }
    lastLayer = layer;
    firstPass = false;
    layer.on({
        click: function(){
            layer.setStyle({
                        radius: 14,
                        weight: 2,
                        opacity: 0,
                        fillColor: '#FFCC00',
                        fillOpacity: 0.8
                    });
        }
    })**
}

**// Reset any highlights after mouseout
function resetHighlight(lastLayer)
{
    layer.resetStyle(lastLayer);
}**


texts: ET Geowizards: Table Already Exists Log Error

texts:
I have a license for ET Geowizards and I keep coming up with an error using various tools. Most recently I have been using the Global Snap Polylines tool. I ran it a few times and then it randomly threw the follwoing error in the log file:

!Start  Snap Polyline Layer Wizard  Date: 2/27/2013  Time: 3:37:44 PM
The table already exists.
ETGW 10.2 for 10.1, 21 July 2012; ArcGIS 10.1.3035:   ET_ExportToPGDBFCDropZM w =4
!Error:  Snap Polyline Layer Wizard  Date: 2/27/2013  Time: 3:37:44 PM

The table does not already exist. I've tried using different output names and even different output workspace locations and the error still persists. I'm wondering if I need to clear the ET Geowizards cache somehow or maybe it's another issue altogether? If anyone knows what is causing this, I would like to hear.

Thanks,
Mike


texts: Scientific approaches for definition of smallest objects

texts: Are there any scientific approaches for defining the smallest recognizable object? For example, a lake has a complex shape, and at a certain scale it is impossible to distinguish its characteristic features and so it doesn't makes sense to display on the map. How to define this scale?


texts: Openlayers 3, getfeatureinfo and geowebcache

texts:
I am having a bit of trouble using the getGetFeatureInfo function on the TileWMS source in OL3. It works, but it represents the X/Y-parameters as I/J in the generated URL.

/geowebcache/service/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=mylayer&LAYERS=mylayer&SRS=EPSG%3A32633&CQL_FILTER=TIME%20TEQUALS2015-12-07T18%3A00%3A00)&INFO_FORMAT=text%2Fhtml&FEATURE_COUNT=50&**I=162&J=125**&WIDTH=256&HEIGHT=256&CRS=EPSG%3A32633&BBOX=294656%2C6706272%2C316320%2C6727936&callback=

This gives me the following error:

  400: The parameters for x and y must both be positive integers.

If I exchange I for X and J for Y in the above URL, it works fine.

How can I make this work? Is the some setting in geowebcache that can be used for it to accept I/J? Or can I do something in OL3 to make it generate X/Y instead of I/J?


texts: No Landsat across northern Australia in 1974-1977

texts:
Why there is no Landsat imagery available for the Gulf of Carpentaria, northern Australia, in 1974-1977?

I have obtained imagery via USGS Earth Explorer from 1972 through to 2017, but there is a gap in 1974-1977, and I cannot seem to find a reason for this. Theoretically, Landsat 1 MSS (1972-1978) and 2 MSS (1975-1982) were around then, but perhaps they just missed the Gulf?


texts: How to resample a raster considering a minimum number of valid pixels

texts:
I'm looking for a function to resample a raster that can consider a minimum number of valid pixels in order to compute the new pixel value.
E.g. If I have to resample a raster of 300m/px to 1Km/px I would consider a window of 3x3 pixels to compute the average value. But I would like to set a control to be sure that at least 5 pixels inside my matrix have valid data.

I'm trying using GDAL (in PyQGIS) but neither GDAL Translate nor GDAL Warp have this option


texts: Flat polygon feature to serve floors to time bins on Z axis of ArcScene?

texts: How can I create square polygons with different z-values to create the horizontal panes shown in the image from ArcScene?


texts:
I have very simple shapefile, which consists of couple of polygons, and I wanted to snap boundaries automatically. Unfortunately qgis doesn't have such function, but grass has according this answer: How to check Topology in QGIS?

I followed provided walk-through, i.e. I made new grass mapset, then imported shapefile layer from qgis, but when I wanted to apply v.clean.break or v.clean.snap command, I couldn't select input vector map as there was not option to select one (FYI I also did refresh in Browser tab):

Browser shows that there is one, but all values are zero, so I wonder if I made a mistake somewhere:

I then deleted created mapset and closed qgis, then opened grass to do the same task but now by just using grass. I imported shapefile to vector mapset but it's just same issue, it's like mapset is empty or something like that.

I use up to date OSGeo4W environment.

What could be issue here?

Update:

(Tue Dec 11 17:09:11 2012)
v.in.ogr dsn=C:\Temp\geo\1.shp layer=1 output=x
Projection of input dataset and current location appear to match
Layer: 1
Counting polygons for 11 features...
Importing map 11 features...
-----------------------------------------------------
Building topology for vector map <x_tmp>...
Registering primitives...
0 primitives registered
0 vertices registered
Number of nodes: 0
Number of primitives: 0
Number of points: 0
Number of lines: 0
Number of boundaries: 0
Number of centroids: 0
Number of areas: -
Number of isles: -
Building topology for vector map <x>...
Registering primitives...
0 primitives registered
0 vertices registered
Building areas...
0 areas built
0 isles built
Attaching islands...
Attaching centroids...
Number of nodes: 0
Number of primitives: 0
Number of points: 0
Number of lines: 0
Number of boundaries: 0
Number of centroids: 0
Number of areas: 0
Number of isles: 0
(Tue Dec 11 17:09:13 2012) Command finished (1 sec)


texts:
[A couple of other questions deal with this sort of problem (see Comparing projections of two datasets when SRID not available and How to test if two GDAL datasets are in the same projection?) but not from an R perspective.]

I'm trying to figure out the best way to compare projections in R because I have a list of shapefiles and I want to reproject them to WGS84 if they are not currently in WGS84.

I thought of several ways to do this, but none seem to be bulletproof:

The simplest way I can think of to do this is to use the equality operator to directly compare two proj4 strings but this will fail, for example, if there are spurious spaces in the strings or differences in decimal places.
I had also thought about trying to determine the EPSG code of a proj4 string and comparing it to 4326 but (I think) not all proj4 strings have corresponding EPSG codes.
Finally, I thought about just using spTransform() on all my layers (~10,000 of them) regardless of whether they're in WGS84 already or not but I figure this will eat a lot of unnecessary compute-time.

I've probably missed something, but I don't seem to be able to find any functions in packages like sp or rgdal that would do a comparison for me.

So, what have I missed? Are there best practices for performing projection comparisons in R?


texts: Determining flood extents based on a daily flow?

texts: I have a DEM raster and daily average flows (e.g., 50Ml/day). And I want to determine how far the water extent goes when the flow accumulates for 7 days based on 50 ml per day?


texts: Arcpy buff does not provide the output attribute "Shape_Area"

texts:
I made a loop to buffer several times a certain network using arcpy. What I noticed is that when I run the code, the output files do not have the attribute fields "Shape_Area" and "Shape_Length" that I get if I perform it directly via ArcMap with the geoprocessing tool "Buffer". Still, the area is shown on the map (therefore its Shape_Area and Shape_Length cannot be zero). Is there an easy way to solve this problem, since I need the "Shape_Area" parameter to work with? My buffer code is simply:

#Enter the name of the fields I want to use to buffer (they exist in "network")
Buffers = ['Name1','Name1','Name1','Name1','Name1']

env.workspace = "C:\Mydirectory"

# Enter the network to be buffered. It contains the fields listed in "Buffers"
network = "FileinMyDirectory.shp"

out = "C:\Anotherdirectory.gdb"
for buff in Buffers:
       arcpy.Buffer_analysis(network, os.path.join(out,buff), buff,)


texts: Which layer for least cost path analysis

texts: New QGIS user here, so maybe it's kind of a stupid question. I want to calcualte the least cost path between two points but I don't know which layer to use. I have two layers, one is a hillshade layer and the other one a color relief map. I tried both and both worked, but which one si the correct layer to use?


texts: How to not allow drag map out of window with OpenLayers 3?

texts:
I can drag the map out of window.

See this example. Can I disable these borders not allowing drag the map out of window?


texts: Eliminate function does not eliminate required boundary/arc?

texts:
Eliminate can eliminate either the common boundary with the adjacent largest polygon or the longest shared boundary.

I want to eliminate/delete common boundry with the third one besides these two.

Can the Eliminate tool do this or is there another tool I need to use?


texts: Visualising most likely direction of migration given points with ages using R?

texts:
I have a question about how best to visualise the predicted direction of migration given a set of points with estimated ages.

Imagine I have a set of geo coordinates, representing communities.  For each point, I have an estimate of the age of that community.  I want to draw an arrow on the map that indicates the most likely direction of general migration over time.  I also want the magnitude of the arrow to represent how consistent the gradient is / how confident we are that there's a meaningful alignment between time and space.  I'm not asking about the theory of the process of migration, just about how you would decide on the angle, position and magnitude of the arrow.

For example, below red points are older than green points.  In the first box, it looks like a sensible conclusion that migration was South-East.  In the second box, there's no clear pattern, so the arrow is smaller.

Is there a standard, principled way of doing this?  One way I thought of was to draw a vector between the oldest point and the 2nd oldest, then between the 2nd oldest and 3rd oldest, and so on.  Then sum the vectors (but how to decide the start location?).  Or maybe you need to take the vectors between all pairs of points and weight them by the relative ages?  Or maybe I need to work out the contours, then draw a path that starts at the highest point and heads downhill?

An R-based solution would be ideal.


texts: Extracting the coverage of reduced resolution

texts:
I have the following reprojected dataset:

//The reprojection variables
var crs = 'EPSG:3857';
var scale = 30000;

//import Datasets
var MOD = ee.ImageCollection('MODIS/006/MOD10A1')
.select('NDSI_Snow_Cover').filterDate('2017-07-01')
.map(function(img){ return img.gte(40).clip(C_A)});
Map.addLayer(MOD,{},'MODIS_FUll')

var MOD_rp = MOD.map(function(img){var rr = img.reduceResolution({
  reducer:ee.Reducer.sum().unweighted(),maxPixels:15000,
bestEffort:false}).reproject({crs:crs,scale:scale});
  return ee.Image(rr)});

I would like to mask all pixels with a coverage less then 30%. This information is visible when I click pixels in the inspector tab:

However, I cannot find it in the properties. Is there a way to extract this info?


texts: PostGIS query points based on distance from regular grid

texts:
I have a PostGIS table with a large number of points (>150k) on a roughly regular grid -- all points are all about 5km apart. Each point has a geometry field and an associated value nr_days.

Is there a clever way to query a subset from this table on a coarser grid, say only points that are 10, 15, 20...km apart?

So far I have only found cases where people were interested in features within a certain distance rather than the other way.


texts: Computing mean of all rasters in a directory using python

texts:
I don't have ArcGIS (ArcPy) so I am hoping to use some other python package such as numpy to compute the mean of a list of rasters in a directory (i.e. Raster1.tif + Raster2.tif + Raster3.tif / 3).

I like to input the directory in the console for computations- using this:

import os
import sys

directory = sys.argv[1]

rasters []
for i in os.listdir(directory):
    if i.endswith(".tif"):

After this I am not really sure what package or code to use. The only examples I have found online are using ArcPy.

Can anyone help me out?

I am hoping to export a single averaged raster as a 16-bit tif.