texts:
I'm working with python and I have to download the whole year of Sentinel-3 data. I'm using sentinelsat library for download data. I have a question about downloading Long Term Archive (LTA) Copernicus products.
Can I make a request for a download for a whole year or can I only make requests for 1 product per 30 minutes?


texts: How do I implement a thematic category theme for Mapsui?

texts:
I want to implement a theme based on category values for a field within an Arcview shapefile.  There is a something similar in SharpMap but not in Mapsui.

Any ideas?


texts:
Possible Duplicate:
  How to import DWG in QGIS project?

I have a .dwg file and I need to covert to shp file using QGIS. There is the dxf to shp plugin, but I need .dwg file to shp file.

How can I do this?


texts: Checking for geospatial similarity in Python?

texts:
Say I have two pandas dataframes in a format similar to this one where each 'Point' represents a person and where they were at. Assume that there are multiple line strings representing different people walking in a specific path. I want to compare based on foot traffic between the two days, whether there is 'similarity' (e.g., based on deviation parameter that could be set) between those of one day vs another one.

Name, Latitude, Longitude, Date
Jordan,<lat>,<lon>, 2017-08-01T00:00:05
Jordan,<lat>,<lon>, 2017-08-01T00:00:08
Jordan,<lat>,<lon>, 2017-08-01T00:00:10
Jordan,<lat>,<lon>, 2017-08-01T00:00:16
Sarah,<lat>,<lon>, 2017-08-01T00:00:20
Sarah,<lat>,<lon>, 2017-08-01T00:00:30
Jordan,<lat>,<lon>, 2017-08-01T00:00:32

I use shapely to construct paths/lines that represent where each person was at on a given day and time.

How I generate the lines...

dayonegeom = [Point(ab) for ab in zip(dayonedataframe.longitude, daytwodataframe.latitude)]

dayonegeodataframe = GeoDataFrame(dayonedataframe, geometry=dayonegeom)
daytwogeodataframe = GeoDataFrame(dayonetwoframe, geometry=daytwogeom)

What is the best way for me to filter the dataframe or GeoDataFrame such that only the paths that are the most 'similar' to each other are kept, while eliminating the ones that are not?

Looking for the best way to do this, be it in pandas before the data gets converted to a geodataframe, or geopandas after it's been converted.


texts: Query Data in ArcOnline REST API. Set-up error

texts:
I am attempting to query data that has been loaded into ArcGIS in the form of layers.

I am attempting to query out basic data (Any data at this point) from a imported Layer.  I think the error is with my .Net code and how I set up to attempt the query but I seem to run into errors when I attempt to query data that is present in the tutorial.

Code Below:

        public async static Task QuerySantaAnaData()
    {
        string tokenString = ""; //This is generated on the the Dashboard for ArcGIS for Developers page under Temporay Access Token

        string path = @"https://services3.arcgis.com/{My ID}/ArcGIS/rest/services/{Service Name}/FeatureServer/0/query?where=FID>0&outFields=*";
        try
        {
            using (HttpClient client = new HttpClient())
            {
                client.Timeout = new TimeSpan(0, 1, 0, 0);
                client.DefaultRequestHeaders.Add("Token", tokenString);

                using (HttpResponseMessage response = await client.GetAsync(path))
                {
                    using (HttpContent content = response.Content)
                    {
                        string data = await content.ReadAsStringAsync();

                    }
                }
            }
        }
        catch (Exception e)
        {
            Console.WriteLine(e.Message);
        }

        Console.ReadKey();
    }

The Inner Exception I receive is

{System.Net.WebException: The underlying connection was closed: An unexpected error occurred on a send. ---> System.IO.IOException: Authentication failed because the remote party has closed the transport stream.
   at System.Net.TlsStream.EndWrite(IAsyncResult asyncResult)
   at System.Net.PooledStream.EndWrite(IAsyncResult asyncResult)
   at System.Net.ConnectStream.WriteHeadersCallback(IAsyncResult ar)
   --- End of inner exception stack trace ---
   at System.Net.HttpWebRequest.EndGetResponse(IAsyncResult asyncResult)
   at System.Net.Http.HttpClientHandler.GetResponseCallback(IAsyncResult ar)}

I have even attempted to perform this API request with the example data that ArcGIS provides in their tutorials here https://developers.arcgis.com/labs/rest/query-a-feature-layer/


texts:
I have found XYZ layers seem to ignore projections or can't seem to use them. I want to operate using a projection so I can buffer point data, for example. I was attempting to use ESRI:102171 because my data is from Victoria. But even using EPSG:3857 and setting as a projection re-projects my other data but does not effect XYZ layers. In other words the basemaps are in one spot and my data is elsewhere. If I return my layers to CRS EPSG:3857 they all line up. Can XYZ tiles not work with projections?

The problem started when, after I had the map set up perfectly using EPSG:4283 (GDA94 - Australia), I tried to buffer my point data by 50m. As the project was in a CRS and not projection the only option was to use degrees when buffering. I then used ESRI:102171 projection to resolve this. XYZ layers didn't come to the party, so to speak. So, I can have a nice map with XYZ tiles and not be able to perform buffer and other operations OR perform these operations in meters without XYZ tiles to support my work. This is my question, can I have both?

I also tried working from EPSG:3857 and using this for projections which had my data lining up but XYZ tiles not complying. QGIS automatically projects on the fly as far as I understand. I just can't get data to match up with XYZ tiles once projected into a local projection.


texts: How to calculate nearest line from point using PostGIS?

texts:
I have point data, which store the "pharmacy_name" field

I also have street data (as linestring), which store the "street_name" field.

I want to know which pharmacy is located on which street. How can I calculate this with a SQL query in PostGIS?


texts: Adding list of values in ModelBuilder using ArcGIS Pro

texts:
I'm searching for the option to use lists in a model of the ModelBuilder under ArcGIS Pro. In ArcMap 10.X, I could right-click on the variable - properties - and under "this variable contains" change the toggle to "A list of values". I cannot find a similar function in ArcGIS Pro.
I explicitly do not search for iterators!
In ArcMap it looked like this:


texts: Leaflet/MarkerCluster plugin - how do update the HTML content of a DivIcon?

texts: In Leaflet, I am using the markercluster plugin, in which you can use the DivIcon class to display custom icons for the clusters. You can set your own HTML for these clusters. However, I can't find how to update the HTML of an existing cluster icon. How can I do this?


texts: Is it possible to disable editing for certain layers in a QGIS project?

texts:
I want to disable editing some (but not all) layers in one QGIS project, while still having them editable in a different project.

Is it possible to disable editing certain layers within a QGIS project without making the actual data read-only?


texts: Uncompressed vs Subdivided: which one offers better PostGIS INTERSECTS performance?

texts:
In my PostGIS I have a point table (roughly 51 millions rows) and a multipolygon table. I need to find a polygon label/ID for every single point in my point table (roughly 51 millions rows) in PostGIS.  So this is considered as point-in-polygon query (with ST_INTERSECTS).

This article suggests that using external storage for geometries significantly improves point-in-polygon queries. Here, geometries in external storage are not compressed hence the performance improvement; but there's potential extra storage size penalty (which is acceptable).

Another way to improve such query is using ST_SUBDIVIDE as explained here. This way, the polygons are divided into rectangular pieces. Point-in-polygon queries, both hit and miss, are improved.

Question:

Uncompressed vs Subdivided, generally which one is more preferable in terms of performance?
If this is rather case-specific, what are the considerations?

Yes I can implement both in a polygon (or multipolygon) table. But when I choose to use small number for max_vertices (a param in ST_SUBDIVIDE), the geometries won't be big enough to bleed to be uncompressed.

Additional infos

polygon complexity:

count    =   319 rows
min      =   97 vertex
median   =   3,942 vertex
max      =   217,436 vertex
stdev    =   22,921.20793 vertex
average  =   10,458.48903 vertex
37 out of 319 have Npoints < 512, so these won't be toasted.

So i made a test using 4 variations for the polygon table:

Conventional
Uncompressed (but not subdivided)
Subdivided (with max vertice= 64, the 319 explodes to 109,574 rows )
Subdivided (with max vertice= 512, the 319 explodes to 19,512 rows )

I tested EXPLAIN ANALYZE with a portion (±500k rows) of my point table (total 51 millions of rows) against those 4 variations. This is to get the 'Execution time'.

EXPLAIN ANALYZE
WITH poi AS (
SELECT  id, geom
FROM    point
WHERE   kpr_id = '430'
)
SELECT  p.id, b.id_polygon
FROM    poi p, polygon b
WHERE   ST_INTERSECTS(b.geom, p.geom) --both geom are indexed
        IS TRUE;

Interestingly, their Execution time are:

Conventional: 9 minutes
Uncompressed (but not subdivided) : 10 minutes
Subdivided (with max vertice= 64) : 8.9 hours (!)
Subdivided (with max vertice= 512) : 1.7 hours

So it seems that using external (to get uncompressed geometry) does not improve performance. ST_Subdivide even gets worse execution time.

Below is the Explain analyze result (with polygon variation #3, subdivided):

Nested Loop  (cost=467106.10..15291893745.32 rows=18183694 width=15) (actual time=8217729.613..32306920.503 rows=526095 loops=1)
  Join Filter: (((b.geom && p.geom) AND _st_intersects(b.geom, p.geom)) IS TRUE)
  Rows Removed by Join Filter: 57645807435
  CTE poi
    ->  Bitmap Heap Scan on point  (cost=5958.78..467106.10 rows=497847 width=36) (actual time=587.892..14814.902 rows=526095 loops=1)
          Recheck Cond: ((kpr_id)::text = '430'::text)
          Rows Removed by Index Recheck: 1942446
          Heap Blocks: exact=105737 lossy=67782
          ->  Bitmap Index Scan on point_kpr_id_idx  (cost=0.00..5834.32 rows=497847 width=0) (actual time=547.728..547.728 rows=526095 loops=1)
                Index Cond: ((kpr_id)::text = '4303130000'::text)
  ->  Seq Scan on polygon b  (cost=0.00..11434.74 rows=109574 width=594) (actual time=5.819..1415.867 rows=109574 loops=1)
  ->  CTE Scan on poi p  (cost=0.00..9956.94 rows=497847 width=36) (actual time=0.016..97.227 rows=526095 loops=109574)
Planning time: 0.415 ms
Execution time: 32306984.195 ms --this is 8,9 hours(!)

And below is all Execution time for all 4 polygon variations with each CTE and non-CTE queries.


texts: Custom Form with QTDesigner: Combine signals and slots?

texts:
I dare to ask that question here and not in the programmers section, cause i'm actually not looking for a coding answer, but a solution through the QtDesigner GUI (if possible).

Custom QT-forms combined with QGIS can be a really powerful solution for data-entry tasks, so i try to tweak things further, but already reached my limits.

Little explanation on the form:

I have lists of bird observations with start and end times.
I already connected some time fields through signals and slots, so when changing the start field, the end fields take over that value as well (only have to change the minutes afterwards -> less typing -> less errors)

When an entry is finished and the next item is digitized, the previous values will be loaded as default.

As most of the time intervals are contigous, i would like to be able to
copy the former end-value to the new start-value when pressing the button E -> A.

In the example that means, that the reloaded end-value 15:03 would be inserted instead of 15:00.

It's quite easy to send values form here to there via signals and slots, but i didn't understand so far if it's possible to trigger such an action also with a third partner (the button)


texts: Using an UpdateCursor inside of a For loop

texts:
I am trying to take PLSS data from a feature class and export it to a Word document. Part of that is cycling through through the feature class and concatenating several text fields into one. Right now, the concatenation works fine, but updating the fields via an UpdateCursor is causing me problems. As it stands, it writes the info from the first record into all the fields instead of cycling through. I see that the two variables legString and infString get set before the UpdateCursor--how can I adjust my code so that I can update those two fields as it goes?
import arcpy
from functools import partial

# Set variables
fc = **PATH TO MY GDB**  # I selected the PLSS_SecondDiv that intersected a project area and exported them
newFields = ['legal_text', 'info_text']

"""Since all arguments are the same except field_name, create partial that has default arguments for everything
except field_name"""
addField = partial(arcpy.AddField_management, fc, field_type="TEXT")  # This lets you create multiple fields in one line

# Loop through list of fields, calling partial to add fields to table
for field in newFields:  # For each field in the field list
    addField(field)  # Adds the field to the table, using Partial()
    arcpy.AddMessage("{} created.".format(field))

# Parse quarters--goes through the feature class, selects the 3 fields, and stores them in "lines"
lines = [row for row in arcpy.da.SearchCursor(fc, ("PLSSID", "FRSTDIVID", "SECDIVLAB"))]  # fc is PLSS feature class
card = ["NE", "NW", "SE", "SW"]  # Cardinal directions list; helps distinguish Lots from Quarter sections

print("Parsing quarters.")
for row in lines:  # for each item in the SearchCursor
    p = row[0]  # PLSSID; used to parse township/range
    f = row[1]  # FRSTDIVID; used to parse section
    s = row[2]  # SECDIVLAB; used to parse lot or quarter-section
    first = s[0:2]  # first 2 characters of SECDIVLAB
    second = s[2:len(s)]  # last 2 characters of SECDIVLAB
    twn = p[5:7]  # 6th and 7th character of PLSSID "21"
    twndir = p[8]  # 9th character of PLSSID "N"
    rng = p[10:12]  # 11th to 13th character of PLSSID "05"
    rngdir = p[13]  # 14th character of PLSSID "E"
    sec = f[-3:-1]  # 4th character from the end to 2nd character from the end of FRSTDIVID

    # Concatenate strings: legString is for legal description, infString is for info block
    if first in card:  # if the first 2 characters are in the card list (means it is a quarter-section)
        legString = '{} ¼ of the {} ¼, Section {}, Township {}{}, Range:{}{}'.format(first, second, sec, twn, twndir, rng, rngdir)
    else:  # if the 2 characters aren't in the card list (means it is a Lot)
        legString = 'Lot {}, Section {}, Township {}{}, Range:{}{}'.format(second,sec,twn,twndir,rng,rngdir)
    infString = 'T/R/S: {}{}/{}{}/{}'.format(twn,twndir,rng,rngdir,sec)  # info block text
    print(legString)
    print(infString)

    #  UpdateCursor to write the strings to fields
    with arcpy.da.UpdateCursor(fc, newFields) as uCursor:  # This is the problem. The two variables are set and this loops
        for cur in uCursor:
            cur[0] = legString
            cur[1] = infString
            uCursor.updateRow(cur)  # No error, but doesn't loop either. Writes 1st record's values to all of them```


texts: Receiving NULL values from Point Sampling Tool in QGIS

texts:
I am using QGIS 2.8.2. My project includes a raster layer of flow depth representing the output from a flood simulation (.DAT file). It also includes a vector layer of points representing the centroids of buildings within the inundated area. I am simply trying to query the raster pixel values (flow depth) at the building centroid locations using Point Sampling Tool, but am receiving NULL values for every entry.

I cannot attach the project files, as it is private data, but a screengrab is here:

The CRS is the same for both the raster layer and point layer (EPSG:2452. This is also set as the CRS for the project, seen bottom-right of the raster layer). The CRSs were set and checked using the "Coordinate Reference System Selector" form accessed by right-clicking the layers and selecting "Set Layer CRS".
The two layers definitely overlap (see screengrab, the black points show the building centroids).
The raster layer is correctly selected using the point selection tool form (see screen grab. There is only one raster layer and one point layer, so no possibility of error).
The raster .DAT file imports fine (using Layer -> Add Layer -> Add Raster Layer), and visualizes correctly (pixel values are correct, band rendering displays the pixel values correctly-see screengrab-, the properties histogram for the layer shows the correct distribution of pixel values).
The output layer from the Point Sampling Tool has the correct number of entries (19,814) so the layer containing the sampling points (building centroids) appears to be correctly used.

Previous queries regarding this problem have all been solved by checking that the CRS is the same for both layers. They are the same in my case, and so my question is:

If the CRS is the same for the sampling point layer and raster layer, and the layers are shown to overlap, why is the point sampling tool returning NULL values for every point?


texts: Meaning of Spatial resolution (m) = 2 rng?

texts:
I'm reading a catalog about SAR sensors.

On page 17, it is said that:

I've read the PolSAR data format for UAVSAR before, But not the other documentations related to InSAR, because I didn't need to.
And I'm really wondering

What does rng mean?
Is it an acronym for a phrase? Then which?


texts: How to load the geometry field of a feature class to PostgreSQL?

texts:
I am trying to load a polygon feature class in a file geodatabse to a PostgreSQL table (extension PostGIS installed). This the process I'm following:

I create table in PostgreSQL with required fields
Read feature class using GeoPandas
Loop through geodataframe containing the polygons from the feature class and insert into postgres table.

I seem to be having problems when inserting the geometry field.

This is the code to create the PostgreSQL table:

import os, csv, logging, sys, time
import psycopg2
import geopandas as gpd

#this is the path to the feature class
fc_vml_polygs = .......\geodatabase.gdb\test'

conn_string = "host=XXXXX" + " port=" + "5432" + " dbname=postgres" + " user=postgres" + " password=postgres"
con = psycopg2.connect(conn_string)
cursor = con.cursor()

create_table = """CREATE TABLE public.vml_polygs
(
    objectid integer NOT NULL,
    fcode smallint,
    shape_length double precision,
    shape_area double precision,
    geom geometry(MultiPolygon,27700),
    CONSTRAINT vml_polygs_pkey PRIMARY KEY (objectid)
)
WITH (
    OIDS = FALSE
)
TABLESPACE pg_default;

ALTER TABLE public.vml_polygs
    OWNER to postgres;

CREATE INDEX vml_polygs_geom_geom_idx
    ON public.vml_polygs USING gist
    (geom)
    TABLESPACE pg_default;"""

cursor.execute(create_table)
con.commit()

Now I read feature class in geodatabase with geopandas and try to load to postgres

#Read polygons feature class
vml_polygs = gpd.read_file(r'........\geodatabase.gdb', driver='FileGDB', layer='test')

#Loop through polygons to load them into pstgres
for index, row in vml_polygs.iterrows():
    load_to_pg = "INSERT INTO public.vml_polygs \
                  VALUES ({0}, {1}, {2}, {3});".format(index, row[0], row[1], row[2])
    cursor.execute(load_to_pg)
    con.commit()

It throws the following error:

Traceback (most recent call last):
  File "footprint.py", line 70, in <module>
    cursor.execute(load_to_pg)
psycopg2.ProgrammingError: syntax error at or near "306481.5500000007"
LINE 1: ...816, 76.65149999882934, MULTIPOLYGON (((452640.71 306481.550...

If I don't try to load the geometry field it works fine. It seems I need to do soemthing with the geometry field of the geodataframe before loading it into PostgreSQL.


texts:
I want a global file chopped into multiple netcdf files corresponding to the modis tiles.
I currently have a directory of GLOBAL netcdf files over time.
NOTE: ERA interim data is on the regular lat/lon grid.

Longitudes range from 0 to 360, which is equivalent to -180 to +180 in Geographic coordinate systems. LINK

This is how my directory currently looks:
ERAIN_SFC00_6H_10mWind.nc  ERAIN_SFC00_6H_2mtemperature.nc  ERAIN_SFC00_6H_TCW-O3.nc

What I want is a set of files on the sinusoidal projection, one file for each MODIS tile. Like this:
ERAIN_SFC00_6H_TCW-O3_h07v06.nc  ERAIN_SFC00_6H_TCW-O3_h24v03.nc         ERAIN_SFC00_6H_10mWind_H15V05.nc  ERAIN_SFC00_6H_10mWind_H31V06.nc
ERAIN_SFC00_6H_TCW-O3_h08v03.nc  ERAIN_SFC00_6H_TCW-O3_h24v04.nc         ERAIN_SFC00_6H_10mWind_H15V07.nc  ERAIN_SFC00_6H_10mWind_H31V07.nc
ERAIN_SFC00_6H_TCW-O3_h08v04.nc  ERAIN_SFC00_6H_TCW-O3_h24v05.nc         ERAIN_SFC00_6H_10mWind_H16V05.nc  ERAIN_SFC00_6H_10mWind_H31V08.nc
ERAIN_SFC00_6H_TCW-O3_h08v05.nc  ERAIN_SFC00_6H_TCW-O3_h24v06.nc         ERAIN_SFC00_6H_10mWind_H16V06.nc  ERAIN_SFC00_6H_10mWind_H31V09.nc
ERAIN_SFC00_6H_TCW-O3_h08v06.nc  ERAIN_SFC00_6H_TCW-O3_h24v07.nc         ERAIN_SFC00_6H_10mWind_H16V07.nc  ERAIN_SFC00_6H_10mWind_H31V10.nc
ERAIN_SFC00_6H_TCW-O3_h08v07.nc  ERAIN_SFC00_6H_TCW-O3_h25v02.nc         ERAIN_SFC00_6H_10mWind_H16V08.nc  ERAIN_SFC00_6H_10mWind_H31V11.nc
ERAIN_SFC00_6H_TCW-O3_h09v02.nc  ERAIN_SFC00_6H_TCW-O3_h25v03.nc         ERAIN_SFC00_6H_10mWind_H17V02.nc  ERAIN_SFC00_6H_10mWind_H31V12.nc

Pseudo code
1. Convert the files to sinusoidal projection
2. determine the extent of the MODIS tiles
3. chop the MODIS tile extents and save as new netCDF

I have absolutely no idea of a) how to convert to sinusoidal projection b) how to chop the modis tiles.
I made these tile extents from already existing modis netcdf files. I don't know if they're accurate (is there an official source defining MODIS tiles?) but they have values.
The data
EXAMPLE DATA HERE
import xarray as xr

ds = xr.open_dataset('ERA_test.nc')
ds

Out[]:
<xarray.Dataset>
Dimensions:        (lat: 256, lon: 512, time: 5)
Coordinates:
  * lon            (lon) float64 0.0 0.7031 1.406 ... -2.109 -1.406 -0.7031
  * lat            (lat) float64 89.46 88.77 88.07 ... -88.07 -88.77 -89.46
  * time           (time) datetime64[ns] 2008-01-02 ... 2008-01-03
Data variables:
    temperature2m  (time, lat, lon) float32 ...
Attributes:
    content:  ERA INTERIM data of var167; converted from grb to netCDF with C...
    history:

 The tile extents
import pickle

tiles = pickle.load(open('tile_extents.pkl','rb'))
tiles['h16v06']

Out[]:
{'lon_min': -23.09401076758503,
 'lon_max': -10.641777724759121,
 'lat_min': 20.0,
 'lat_max': 30.0}

Converting to sinusoidal projection
This is the one that I really don't know how to do. I think I have found the MODIS projection from pyproj.
import pyproj as Proj

p_modis_grid = Proj('+proj=sinu +R=6371007.181 +nadgrids=@null +wktext')

But how do I go about using this to reproject the netCDF data?
Happy to know if I'm using the wrong tool.
Will rasterio be better suited?
Or NCO/CDO?


texts:
I am trying to overwrite an existing fgdb table with a R-script and the arcgisbinding package. I want to then place this script in Model Builder and use the resulting table in the rest of the model.

Here is a simplified version of the R-Script.

tool_exec <- function(in_params, out_params) {
  inTable <- in_params[[1]]
  dfOpen <- arc.open(inTable)
  df <- arc.select(dfOpen)
  df$A <- 0
  df$B <- "xxx"
  # Create the output table
  arc.write(inTable, df)

}

It gives the error:

  Error in .call_proxy : The table already exists.

I have tried giving the output table a different name, but then if I run the script again I end up with the same problem.
If I change the Type of the OutputTable parameter from Derived to Required then I can specify a new table name each time (which I dont want to do) but then once I run the tool once, I have to specify the OutputTable location again manually.

Is there a way around this error?

Notes The box is ticked for the Geoprocessing Options - Overwrite the outputs of geoprocessing operations. I have very limited Python knowledge so that is why I am trying to do it in R & Model Builder.


texts: Get username on OAuth Earth Engine authentication

texts:
I'm using the following code to authenticate users in a web application based on Google Earth Engine:

var onImmediateFailed = function () {
    $('.g-sign-in').removeClass('hidden');
    $('.output').text('Please log in to use the app.');
    $('.g-sign-in .button').click(function () {
        ee.data.authenticateViaPopup(function () {
            // If the login succeeds, hide the login button and run the analysis.
            $('.g-sign-in').addClass('hidden');
            loadMap();
        });
    });
};

// Attempt to authenticate using existing credentials.
ee.data.authenticate(CLIENT_ID, loadMap, null, null, onImmediateFailed);

It works perfectly but I'm not able to get the username of the authenticated user.

How can I get it?


texts:
I'm using ArcMap 10.3.1 from a basic desktop concurrent use license. The problem I'm experiencing occurs when I load an existing mxd (creating using 10.3.1), and then try to switch from data view to layout view. When I try to switch views, the program crashes. I have tried to delete the normal.mxt file and this has not solved the problem. I've also updated my video driver and run the ERSI diagnostic to see if my computer is capable of running the program. I've also run a diagnostic to ensure that my computer is not experiencing any issue with memory or drivers.
Any ideas on what might be causing the problem or how to fix it?


texts:
I am trying to setup a custom basemap using a GTiff in the IAU's Moon2000 projection (30100 or 30101).  The basemap image was originally created through gdal_merge and gdalinfo returns the expected, human-readable WKT CRS information.

As per the geoserver documentation, I then used gdal_retile.py to automatically generate my pyramids and the accompanying directory structure.  This was then moved into /data-dir/data.  A new workspace, for organization, and new store were created.  I am able to add the basemap "layer", using the pyramid plugin, but the CRS defaulted a UTM zone (18N).

As per the geoserver documentation, I then edited the epgs.properties file and added the single line WKT from the link above.  Through the geoserver web administration tool, I can see that the CRS was successfully added, but the wkt is empty.  See the screenshot below.

I suspect that I now need to go in an manually add both a custom datum and spheroid.  Is this the case?  If so does anyone have documentation of this process?  Has anyone gotten planetary data (any planetary body) added to geoserver as a basemap?

The end goal is to server the tiled basemap as a WMS to openlayers.


texts: Drawing circle in Google Earth around point?

texts:
I can't seem to figure out how to draw a circle in Google Earth.

I'm trying to draw a 750km circle around Vancouver.

I'm trying to generate the kml file here:
    http://kml4earth.appspot.com/circlegen.html

But I can't seem to figure out the correct corrodents.
I'm using:
    http://www.tageo.com/index-e-ca-v-02-d-m812033.htm


texts:
i have the following problem: i have downloaded tiff files with gMapMaker and transformed them to GeoTiff using the associated .map files in QGIS. now i want to render these GeoTiff files to use them as a map (like using Google Maps API) and then apply some geographic info (as an SVG element on a particular geographic point) for example using d3.js or OpenLayers.. The important thing is that i have to do all of this offline (yes, i'm trying to re-create a sort of Google Maps offline)..
until now, i tried to use Javascript but i didn't find any library to load GeoTiffs; moreover i can use also other languages (i have no restrictions from this point of view)... can you help me?


texts: Using JavaScript to manipulate the display of GeoJSON on Openlayer4

texts:
I am loading a number of geojson data (all polygons) and using 3 sliders, the user should be able to get a polygon displayed based on the slide value. Here is an online example of what I intend to achieve. In the online example, as I move the slider, each risk layer displayed changes in coverage (one increases and the other decreases). I have created an example here with two geojson data (small and agriculture). How do I get to target each of the small.geojson polygon to show or hide based on the slider values.

I can't seem to figure out how to go about this as each small polygons (as in the example link I created shows that the data is now a canvas.

PS: This is more of a guide need.


texts: Are with WebMercator the X and Y extents identical?

texts:
I am calculating the tile coordinates in a WebMercator projection at different zoom levels and have the following problem:

Cartopy suggests that the projection extent in meters is slightly different in the X and Y direction:

In [67]:     extent_x = ccrs.GOOGLE_MERCATOR.x_limits[1] - ccrs.GOOGLE_MERCATOR.x_limits[0]
   ....:     extent_y = ccrs.GOOGLE_MERCATOR.y_limits[1] - ccrs.GOOGLE_MERCATOR.y_limits[0]
   ....:

In [68]: extent_x
Out[68]: 40075016.685578488

In [69]: extent_y
Out[69]: 40075016.685578495

In [70]: extent_x - extent_y
Out[70]: -7.4505805969238281e-09

While this reference (see the __init__ and MetersToPixels functions in the script at the end) suggests that extent is equal in both directions:

self.initialResolution = 2 * math.pi * 6378137 / self.tileSize
...
res = self.Resolution( zoom )
px = (mx + self.originShift) / res
py = (my + self.originShift) / res

Which is right? Is the discrepancy in X and Y extent that results from using the cartopy definition merely a floating-point approximation issue or real?

If it is an floating point artifact, is there a way to obtain correct X and Y extents from cartopy?


texts: QGIS 2.18: Speed up QgsFeatureRequest in Python Script

texts:
I have two input layers. The first one is a line-layer named "graph" and the second one is a point-layer named "nodes".

The "graph"-layer:

The graph layer holds 177350 line features which are representing a streetmap.
Loaded into QGIS it looks like this:

Each street (or line) holds some information: The id of the street (edgeid), the node id of the startnode, the node id of the endnode and some other attributes as shown below:

The "nodes"-layer:

The nodes layer holds 87667 point features which are representing the start- or endnodes of the streets in the graph-layer. Loaded into QGIS and overlayed with the graph-layer it looks like this:

Every node in the point-layer has a unique nodeid, x-, y- coordinates and a type (the type attribute is not relevant):

What I want to do is to calculate the azimuth for every street (line) in the graph layer. The outputlayer of my python script should look exactly like the graph-layer, just with an additional azimuth attribute. Like this:

To achieve this, I have written an python script which loops through all the features of the graph layer. For every street, I look up the x-,y-coordinates of the startnode, and of the endnode. I will need the x-,y-coordinates for calculate the azimuth later on in the script. I am doing this using QgsFeatureRequest. But as there are 177350 lines to go through, I will have to request 177350 * 2 times (177350 times for startnodes and 177350 times for endnodes). And it just takes too long.

I have already tried to set the "NoGeometry"-Flag and the setSubsetOfAttributes-method as it was suggested here. But it seems to be useless.
What am I doing wrong here? How can I speed up the process of searching for the nodeids in the nodes-layer?

i = 0
n = len(feats)

# loop through graph layer and find node coordinates for from/to nodes
for i, feat1 in enumerate(feats):
    progress.setPercentage(int(100 * i / n))
    fromNodeId = feat1.attributes()[idx_from]
    toNodeId = feat1.attributes()[idx_to]

    # search for fromNode by fromNodeId
    exp = QgsExpression('nodeid = \'' + fromNodeId +'\' ')
    request = QgsFeatureRequest(exp)
    request.setFlags(QgsFeatureRequest.NoGeometry)
    request.setSubsetOfAttributes(['nodeid', 'x', 'y'], inputNodes.fields())
    request.setLimit(1)
    print "time 1"

    for feat2 in inputNodes.getFeatures(request):
        from_x = feat2.attributes()[idx_x]
        from_y = feat2.attributes()[idx_y]
        print "time 2"

    # search for specific ToNode by ToNodeId
    exp = QgsExpression('nodeid = \'' + toNodeId +'\' ')
    request = QgsFeatureRequest(exp)
    request.setFlags(QgsFeatureRequest.NoGeometry)
    request.setSubsetOfAttributes(['nodeid', 'x', 'y'], inputNodes.fields())
    request.setLimit(1)
    print "time 3"

    for feat2 in inputNodes.getFeatures(request):
        to_x = feat2.attributes()[idx_x]
        to_y = feat2.attributes()[idx_y]
        print "time 4"

    fromPoint = QgsPoint(from_x, from_y)
    toPoint = QgsPoint(to_x, to_y)

    newLine = QgsFeature(outLayer.pendingFields())
    (...)
    azimuth = fromPoint.azimuth(toPoint)
    newLine.setAttribute(AZIMUTH, azimuth)
    (...)

EDIT: I have added some timing print outputs. From print "time 1" to print "time 2" it takes approximately 4 seconds. The same with print "time 3" to print "time 4". This means one runthrough the outher for loop takes approximately 8 seconds (I know measuring the time with prints isn't the best way. But I just want to show the bottlenecks of the python code).

If you need some more information on my question, please feel free to leave a comment down below!


texts: Load multiple file geodatabase tables into R directly from a FGDB (arcgisbinding package)

texts:
I have a large number of file geodatabase tables that I'd like to load directly into R as data frames. I have upwards of 80 - 100 tables and it would be much more convenient to automate this process.

I've been playing around with the arcgisbinding package in R, specifically the arc.open() and arc.select() functions, and although loading single tables at a time is easy, I'm struggling to write a script that will execute this process for multiple tables at once. Here's my script so far:

library(arcgisbinding)
arc.check_product()

# Assign my fgdb to a variable and create a character vector of the fgdb table names.
fgdb <- "D:/WorkingDirectory/ZStatTables.gdb"
TableList <- arc.open(fgdb)@children$Table

Following code I found here, I tried to loop through the vector of FGDB table names to load each table into R.

a <- NULL
for(i in 1:length(TableList)){
  a[[i]] <- arc.select(arc.open(paste(gdb,"/",TableList[i], sep = "")))
  }

This produces a list of all of the tables, but their names have not been assigned to them. Is there a way to assign each table it's correct name? Or an altogether better way to load a lot of FGDB tables into R?

Any ideas, either using the arcgisbinding package or other methods?


texts:
I are looking to implement a method to fix self-intersecting polygons by editing the vertex order.  We want to take the polygon rings and re-arrange/add in vertices to make a valid polygon.  This is similar to how ArcMap's Repair Geometry runs.  I have used other methods such as PostGIS's ST_MakeValue, but we didn't like the results because it added extra area to the polygon.  It is because of this, that I am looking to implement my own strategy.

I am looking to take the WKT of the feature shape and manually edit the coordinates.  Doing something like a buffer operation, split and union, or anything else won't really work.  What we want to do is be able to detect and add vertices to self-intersections similar to ArcMap's tool.

If I can get an explanation on some popular tools used to correct self-intersecting polygons and perhaps a code sample/algorithm sample, it would help me greatly.


texts:
Related to my Problem described in PyQGIS memory issues, I'd like to try using destructors for QgsGeometry and other classes. I dont get how to use them though.

According to https://qgis.org/api/2.18/classQgsGeometry.html#aacaf2856a136d270dcf274649439adf7 it should be:
geom.~QgsGeometry()

This is invalid syntax though. Are they explicitly callable at all? And if so - how to do it?

(Qgis 2.18)


texts: What is the default unit for slope in QGIS?

texts: When using Raster Based Terrain Analysis on QGIS for creating a slope layer, what units are the used in the output? Degrees or %?


texts: I have an sf polyline layer in R an would like to calculate the sinuosity (i.e. actual length of segment divided by the shortest distance between start and endpoint). Does anyone know of a function from sf or similar packages to accomplish this?


texts:
I created one geodataframe with geopandas. The geometry field was filled with Multipoint. However, when I try to plot the map to see multipoints, my map returned empty.

How can I plot geometry field (multipoint) using geopandas?

Here is a piece of it:

from geopandas import GeoSeries
from shapely.geometry import Point, MultiPoint

gs = GeoSeries([MultiPoint([(-120, 45), (121, 46)]),
                MultiPoint([(-121.2, 46), (-123, 48)]),
                MultiPoint([(-122.9, 47.5), (-125, 49)])])
print(gs)

gs.plot(marker='*', color='red', markersize=12, figsize=(4, 4))


texts: Data lost after clipping raster layer with drawn rectange, not another layer

texts:
I am trying to clip a raster layer which has a CRS of EPSG:4326, WGS 84. The data disappears after I clip it. I am not using a polygon from another layer; I just drew a rectangle to select the area.

Here is the text in the clip window:

gdal_translate -projwin 32.6739583333 15.3614583333 48.2708333333 3.28645833333 -of GTiff "/home/ross/saffron/Climate data/QGIS/New QGIS from WC/wc2.0_30s_tmax/wc2.0_30s_tmax_01.tif" "Jan T max Ethiopia"


texts: QGIS2Web won't display OpenLayers but does display Leaflet

texts:
I'm trying to make a web map. I have everything working fine in QGIS 2.18.15 (all point exist where they should and layers are all good), but when I try to make a web map with OpenLayers nothing shows up, either on the preview or on export. I can get a base map, but none of my shapes or points appear. When I choose the Leaflet option The points and shapes appear, but the points lose their gradation that I've defined in QGIS (more populated points are bigger and have a brighter hue in QGIS, but all points are the same size and colour in OpenLayers.)

Any help?


texts: Copying attributes from large number of raster files in ArcGIS to Excel?

texts: I am new to GIS, I am having large number of raster files(800) in ArcGis and I want to copy the attribute table of each raster file separately to Excel. Can anyone help me how to perform it in batch mode or using coding?


texts: Geocoding table when more rows expected to be appended to table?

texts:
I'm starting on a project with ArcGIS 10 where I'll be doing quite a bit of geocoding addresses.  I understand the basics of how to do that.  But a frustration on an earlier similar project is that when the boss would come down the hall with just a few more addresses for me to add to the map (over and over). So I would have to repeatedly re-geocode the whole table each time I add a few.

Thinking "there must be a better way", I came across the option on the geocoding dialog to create a dynamic feature class related to a table.  This promises to automatically update the feature class when the data in the table changes.  So I have tried this.  I built a file geodatabase for this project, and imported some test data into it, and then geocoded the addresses.  In ArcCatalog, I see my data table, the feature class, and a relationship class that links them.

However, when I edit an address in the table, I don't see any auto-updating going on.  And more importantly, when I append more rows (via the Append toolbox item), they don't get auto-geocoded, either.  Am I missing a step?


texts: I have 5 PostgreSQL_PostGIS database in my company and i conneting on them through ArcSDE. Four of database work correctly, only with one i have i realy big problem. When i make two versions of database, and i edit data in one version, reconcile and post on DEFAULT version, database do some kinde of automaticly reconcile on second version, and i see all the changes in second version before i do reconcile. Does anyone know what could be the problem?


texts: ArCGIS Server 10/10.1 And Sql Server 2008

texts: Is there any compulsion that SQL server 2008 and ArcGIS server 10/10.1 should be installed on same server machine and is it possible to install SQL server 2008 and ArcGis server on separate server machines?If is it possible then what are the complications and performance issue involved in above scenario.?


texts: Does mixing geometries of different granularity affect POSTGIS index efficiency?

texts:
My POSTGIS table contains a lot of small non-overlapping geometries and I do mostly bbox queries. I would like to add larger / coarser grained geometries to the database (counties, provinces, country) and get an efficient bbox result.

Will my GIST-index query performance be affected significantly if I combine geometries of different granularity in the same table? I'm considering to move the larger geometries to a different table and do 2 queries.


texts: If-Then-Else Statements With Rasters in ModelBuilder

texts:
Basically, I have two models that I have built. Their core processes are similar, but they both have different beginning processes. There is one beginning input layer. The cell values within this layer determine which model will need to run. For example, all cells that are < 5 will need to run Model1. All cells >5 will need to run Model2. Almost all of the processes within both sides of the model require Raster Calculator.

I have set-up an If Field Value Logical statement at the beginning, and then divided the two processes. However, It seems that when I do run the model, It only completes the FALSE portion of the model. For the TRUE portion, it states "All the inputs are not current." How can I make it so it runs both? I will have two outputs: Raster1, and Raster2.

I am by no means an expert at ModelBuilder, and only know the very basics of python. So please excuse me if I am missing something obvious.


texts: Using geography cast as geometry in QGIS No SRID in Airports table. Other table have geom and is 26986

texts:
I trying to add a layer into the map in QGIS. It will execute the SQL and show the answer but will not add it to the map. It states there a geometry error. I type cast geom::geography it should allow it.
SELECT name, iso_country, iso_region, ST_Transform(geog::geometry,4326)
FROM ch10.airports
WHERE ST_DWithin(geog, ST_Point(-94.173653,36.067833)::geography, 100000)
LIMIT 5;


texts: Is there a way to snap to apparent intersections in ArcGIS 10.1+?

texts:
Suppose there are two lines in an ArcMap 10.1 or 10.2 editing session that cross but do not intersect, per se. This is the case when new lines a drawn. Is there a way to snap to the point where the lines cross?

AutoCAD has this ability, and calls these intersections apparent intersections. In the situation below, two new lines (Line 1 & Line 2) are drawn. In Example X, Line 2 was just created and topology has not been checked. The cursor is hollow because the create new feature tool is selected.

In Example X, the snapping tip feature indicates there is no node to snap to at Point A. However, in Example Y, there is a point to snap to at Point A, because topology was checked in Example Y.

Example X

Example Y

So, short of checking for topology errors, there does not seem to be a way to snap to apparent intersections.

Conditions notes

All snapping options are on:

Classic snapping is not enabled.
The only relevant topology rules to this layer are:

The system here is running ArcGIS 10.2, with features in 10.1 SDE on Oracle.


texts:
I have been using Leaflet.js but since I wanted the webgl performance I switched to the newly OpenLayers 3 which seems awesome.

I had bounds on my Leaflet.js Map:

var sw = L.latLng(-90, 180),
    ne = L.latLng(90, -180),
    maxBounds = L.latLngBounds(sw, ne);

map = new L.Map('map', {
        attribution: "",
        attributionControl: false,
        maxBounds: maxBounds,
        zoomControl: false,
        zoomAnimation: true
    }
);

where I only had a single projection being allowed to be displayed and if the user would attempt to go out of bounds, the map would bounce back with an animation.

Is there a way to have something similar with OpenLayers 3?

I basically want to prevent dragging or changing center beyond the actual bounds of a single projection.

I used the following code with OpenLayers 3:

map.on('postrender', function () {
    var view = map.getView(),
        worldExtent = view.getProjection().getExtent(),
        extent = view.calculateExtent(map.getSize());

    if (!ol.extent.containsExtent(worldExtent, extent)) {
        var newExtent = ol.extent.getIntersection(extent, restrictedExtent);
        view.fit(newExtent, map.getSize());
    }
});

which displays a signle projection all right but I can drag or move the center out of it.


texts:
I'm trying to create a custom Polygon style for my layer style.

I receive some information about position, orientation and age for a number of objects and I want to represent those as zoom independent (like using Icons or RegularShapes) isosceles triangles, pointing at an angle and color-coded according to age.

I have a good idea of the color-coding and rotation part, but I'm stumped at creating an isosceles triangle using the style objects provided by OpenLayers.

This is what I have so far:

style = [new ol.style.Style({
    geometry: function(feature) {
        return new ol.geom.Polygon([[
            [0, 0],
            [2, 0],
            [1, 3]
        ]]);
    },
    fill: new ol.style.Fill({color: '#000'}) //black right now
})]


texts: Adding related table to ArcGIS Online feature layer

texts:
I'm trying to figure out how to add a relationship class (related table one-to-many) to a feature layer within AGOL. I do not have ArcGIS Desktop or Server (do have QGIS but I don't think that will help here).

The only info I was able to find led me to believe that it may be possible by editing the service definition JSON but I can't find any reference as to how to edit it to make it work.

What are my options here?

Also, is viewing/editing/adding to related tables supported in AGOL web apps yet (not just Collector)?


texts: Get input layer extent as new GeoJSON polygon on the command line

texts:
Starting from a GeoJSON file containing a country boundary polygon, I would like to obtain a rectangular approximation: another GeoJSON consisting of the bounding box/envelope of the original polygon. This seems simple enough to be possible from the command line, so I was thinking of using ogr2ogr (other options also welcome). But I cannot seem to find a command to do this. Can this be done without getting too cumbersome?
With ogrinfo, it's possible to get the layer extent, but it does not generate a geometry file.
The GeoJSON driver also has a WRITE_BBOX option, but that's just additional info and has no impact on the actual geometry.


texts: QGIS 3 crashes when adding a PostGIS view as a map layer?

texts:
I am trying to load a PostGIS view public.view1 into QGIS (3.6) as follows (per Add PostGIS layer to QGIS via Python Console):

tablename = "view1"
geometrycol = "geom"

from qgis.core import QgsVectorLayer, QgsDataSourceUri
uri = QgsDataSourceUri()
uri.setConnection("localhost", "5432", "user", "user", "")
uri.setDataSource ("public", tablename, geometrycol)
vlayer=QgsVectorLayer (uri .uri(False), tablename, "postgres")
QgsProject.instance().addMapLayer(vlayer)

The view is created from a source table:

create view view1 as select * from my_table;

Both the source table my_table and the view view1 can be loaded via the QGIS 3 GUI manually without a problem. But when loading the view programmatically using the above code, QGIS crashes on the last line:

QgsProject.instance().addMapLayer(vlayer)

The same code works fine on the source table. And this is with PostgreSQL 11.3, PostGIS 2.5, and QGIS 3.6.

Has anyone else encountered similar problems and know of a solution?


texts: Drag Vertices - ol2 Modify Control vs. ol3+ Modify Interaction

texts:
In OpenLayers 2, clicking a geometry with an active modify control running, created a drag vertice on every edge of the geometry which could be used to modify the geometry by dragging.

In OpenLayers 3+ with the modify interaction, things are different. Clicking a geometry instantly activates one vertice only that can be dragged.

Can I obtain the ol2 behaviour in ol3+ as well?


texts: Problem with encoding of CSV and SHP files in Linux

texts:
I am new to using QGIS.  I came up against the encoding issue of shp files.
This solution to solved my problem in Windows (for source I had to put windows-1253 instead of utf-8)

My question is what solution would I use in Linux systems where I have the very same issue?  Where can I "insert" a similar choice?


texts: Raster terrain analysis on QGIS 2.0.1 - No curvature analysis?

texts:
I need to perform a curvature analysis on a DEM. I've found this tutorial

http://www.gistutor.com/quantum-gis/20-intermediate-quantum-gis-tutorials/48-quantum-gis-qgis-raster-based-terrain-analysis-techniques.html

that mentions a "Total curvature analysis" on the Raster terrain analysis plugin. However I've been unable to find that option in the plugin's menu. I am relatively new to QGIS, so I may be missing something, or has that function been discontinued?

Thanks in advance


texts: Batch raster calculator with ModelBuilder in ArcMap 10.6

texts:
I have 30 images of Sentinel-2 data.
How can I calculate NDVI in batch?
I tried to solve the problem with ModelBuilder, but there can only be 1 iteration tool in the model. I additionally created a submodel, but each raster from the submodel is iterated over with all from the main one.

And the images don't match. And after first iteration
Float((Raster(r"С:\Input\S2B_MSIL2A_PCL\ PCL _B08_10m.jp2") - Raster(r"C:\Input\S2A_MSIL2A_SXT\ SXT _B04_10m.jp2")) / (Raster(r"С:\Input\S2B_MSIL2A_PCL\ PCL _B08_10m.jp2") + Raster(r"C:\Input\S2A_MSIL2A_SXT\ SXT _B04_10m.jp2")))<br/
I will add the question
<br/Here I'v got 15 subfolders with sentinel data - C:\Input\S2*\S2A*\GRANULE\L2A*\IMG_DATA\R10m\bandX.jp2
I want to understand how it's possible to process all rasters using python.
For example, I put 4 Bands (Band4 x2 and Band8 x2) to C:\Input and use this code
arcpy.env.workspace = r"<C:\Input>" #define the workspace (your input path)

raster_list = arcpy.ListRasters(("*B04_10m"))
for raster in raster_list:
print("YEAH!")
rasRed = Raster(raster)
rasNIR = Raster(raster.replace('B04','B08')) #corresponding NIR band
outraster = float((rasNIR-rasRed)/(rasNIR+rasRed))
outraster.save(raster.replace('B04','NDVI'))#change the input name to NDVI

There is an error in line 4, 'NoneType' object is not iterable.
If I use these two codes - no error and no result.
import glob

 for raster in glob.glob(r"C:\Input\S2*\*B04_10m.jp2"):
     print(raster)
     rasRed = Raster(raster)
     rasNIR = Raster(raster.replace('B04','B08'))
     outRaster = float(rasNIR-rasRed)/float(rasNIR+rasRed)
     outraster.save(raster.replace('B04','NDVI'))

import glob

 for raster in glob.glob(r"C:\Sentinel_input\S2*\S2*\GRANULE\L2*\IMG*\R10*\*B04_10m.jp2"):
     print(raster)
     rasRed = Raster('raster')
     rasNIR = Raster(raster.replace('B04','B08'))
     outRaster = float(rasNIR-rasRed)/float(rasNIR+rasRed)
     outraster.save(raster.replace('B04','NDVI'))


texts: Option to export from ArcMap 10.7 to CAD not in group

texts:
Is there an option to prevent ArcMap 10.7 from exporting CAD polygons as group?
I am using ArcMap to export layers to CAD as DWG. Does someone know how to prevent ArcMap to export polygons as group?


texts: Implementing Trace Geometric Network into Model Builder

texts: I would like to use Trace Geometric Network in model builder with the option to have a model parameter for the flags (like the flags I add from Utility Network Analyst Toolbar) and return the results as a selection set.


texts: How to convert raster to polyline in QGIS?

texts: How can I convert raster (stream segments result of a watershed analysis) to vector polylines, preferably shapefile. In ArcGIS this is easily acomplished with the Raster to Polyline tool.  I only see Polygonize (Raster to Vector) tool under the Raster -> Conversion menu.


texts: For which operations a PostGIS Geography index is useful?

texts:
I need to do KNN- and within-region queries for POI data. I store the POI locations in a geography column. Since I understand that ST_DISTANCE doesn't use the geography index, I added a geometry index in order to use the <-> operator with indexed performance. This works nicely (unfortunately not over the 180° meridian but that's OK in my case).

Now I'm asking myself for which operations a geography index is actually useful?


texts: Vector data(point) is drawn in wrong place on map (OpenLayers)

texts:
I need to show a map (world map, used default OpenLayers WMS) and one point on it (with events like onhover, etc). Here is my code:

                   var options = {
                            projection: new OpenLayers.Projection("EPSG:900913"),
                            maxResolution: 6000
                        };
                        map = new OpenLayers.Map('map', options);

                        var wmsLayer = new OpenLayers.Layer.WMS(
                        "OpenLayers WMS",
                        "http://vmap0.tiles.osgeo.org/wms/vmap0",
                        {layers: 'basic'}
                    );

                    var vectors = new OpenLayers.Layer.Vector("Vector Layer");

                    point = new OpenLayers.Geometry.Point(20.088844299316406, 51.8321709083475);
                    vectors.addFeatures([new OpenLayers.Feature.Vector(point)]);

                    map.addLayers([wmsLayer, vectors]);

                    map.zoomToMaxExtent();

But this code locates the point is drawn not in correct place, but somewhere near Africa (that place lat and lon is 0, 0). Question: Why this happens and can I fix it? I just need to locate the point to the correct place. Paradox when I print this point in console then it shows that the point lan and lon are as needed (as defined). But it is still in the wrong place...


texts: can I export service area (network analysis) to shapefile?

texts: Is it possible to export service area (network analysis) to shapefile? if so, how can I do it? I am using ArcGIS 10.


texts: Removing grey margin from exported map of ArcGIS Desktop?

texts:
In above image there is a gray color margin.

Is there any way I can change the color of it in to white?


texts: Problem creating FeatureLayer object in esri javascript API

texts:
I'm just a simple farmer tending to his ArcGIS online hosted feature layers.  I pinpointed the problem as that I can't seem to properly create a FeatureLayer object based on a url to a hosted feature layer.  I can execute query tasks and what not just fine, but I'm really not sure where I'm going wrong with creating a feature layer.  I tried just using an esri example feature layer, and it's still not working for me.  Am I just missing something dumb?  Javascript API v3.31 btw

My code:

<!DOCTYPE html>
<html>
<head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <title>testfeatures</title>
    <script src="https://js.arcgis.com/3.31/"></script>
    <script>
      require([
        "dojo/_base/declare",
        "dojo/dom",
        "dojo/parser",
        "dojo/ready",
        "esri/layers/FeatureLayer"
    ], function(declare, dom, parser, ready, FeatureLayer) {
        parser.parse();

        ready(function() {
            var testLayer = new FeatureLayer("http://sampleserver6.arcgisonline.com/arcgis/rest/services/USA/MapServer/0", {
                mode: FeatureLayer.MODE_AUTO,
                outFields: ["*"]
            });

            console.log(testLayer);
            console.log(testLayer.fields)
        })
    });
    </script>
</head>
<body>

</body>

Console logging the fields gives "undefined", whereas logging the feature layer itself gives me an object.
Thanks for any help!


texts:
I'm looking for a function/formula/library that basically takes a MGRS square identifier such as 4QFJ16 and returns me the coordinates of the 4 corners of this polygon in some simple form such as EPSG:4326 lat/lon pair.

The idea behind is that I regularly have data to display where we only know in which square it was taken (no northing/easting information). So I'd simply like to be able to display this specific polygon on a map.

What I can find easily is large graticules as shapefiles, or information to convert a specific point in this square to lat/lon coordinates.

But some "as simple as possible" formula would be better, since I could easily make browser calculations instead of embedding a huge amount of data (the complete graticule). Also, a formula would have the advantage to be language-agnostic, so I can use it for different use cases (realtime JS conversion JS, filling a PostGIS database, ...).

I don't need a huge precision, I work with 10 or 1km MGRS squares.

Hope my question is clear. Many thanks in advance.


texts: Drawing WMS layer with Attribute-Based Dynamic SLD to HTML Canvas?

texts:
I'm working on an app using GeoServer and OpenLayers 2.13.1. I have several WMS layers and I need to export whatever in my openlayers map div as an image. After researching, I found this blog: http://blogs.edina.ac.uk/category/openlayers-html5-canvas/ which is working fine (exporting my layers as an image).

However, my wms layer also has dynamic sld functionality (using SLD_BODY). The problem shows up when I use attribute-based sld: it wasn't drawn on the canvas. It only works when I use non-attribute-based sld.

Did I miss something here?

This is how I call the wms layer:

wms = new OpenLayers.Layer.WMS(layer_name,
    WMSServiceURL, {
        layers: "bappeko:"+map_name,
        styles: null,
        transparent: true,
        format: "image/jpeg"
    },{
        ratio: 1,
        singleTile: true,
        tileOptions: {maxGetUrlLength: 2048},
        buffer: 0,
        projection: new OpenLayers.Projection("EPSG:32749")
    }
);
map.addLayer(wms);

This is JS function I used to export:

function exportMap() {
var mapcanvasDiv = null ;
var mapDiv = document.getElementById("mapdiv") ;

if (mapcanvasDiv === null) {
    var canvasElement = document.createElement("canvas");
    canvasElement.id = "mapcvs";
    mapcanvasDiv = document.getElementById("mapcanvas");
    mapcanvasDiv.appendChild(canvasElement) ;

    canvasElement.width = mapDiv.clientWidth ;
    canvasElement.height = mapDiv.clientHeight ;
}

var mapCanvas = document.getElementById("mapcvs");
var mapContainer = document.getElementById("OpenLayers_Map_2_OpenLayers_Container");

if (mapCanvas !== null) {
    var ctx = mapCanvas.getContext("2d");
    ctx.clearRect(0, 0, mapCanvas.width, mapCanvas.height);

    var olLayersDiv = document.getElementsByClassName("olLayerDiv");
    var olLayers = map.layers;

    for (var i=0; i<olLayersDiv.length; i++){
        //ONLY EXPORT VISIBLE LAYERS
        if (olLayers[i].visibility == true) {
            var layertiles = olLayersDiv[i].getElementsByClassName("olTileImage");
            for(var j=0; j<layertiles.length; j++) {
                var tileImg = layertiles[j];
                var offsetLeft = tileImg.offsetLeft;
                var offsetTop = tileImg.offsetTop;
                var left = Number(mapContainer.style.left.slice(0, mapContainer.style.left.indexOf("p")));
                var top = Number(mapContainer.style.top.slice(0, mapContainer.style.top.indexOf("p")));
                ctx.drawImage(tileImg, 0, 0);
            }
        }
    }

    var dataUrl = mapCanvas.toDataURL();
    window.open(dataUrl, "toDataURL() image", "width=600, height=360");
}

This is how I set my dynamic SLD (the SLD is created dynamically by PHP script via Ajax) in JS:

l = map.getLayersByName(layer_name)[0];
l.mergeNewParams({sld_body: sld, styles: null});

I guess I've found the culprit: It seems that my PHP-created SLD is way too long even when I've set tileOptions: {maxGetUrlLength: 2048} because I use many <ogc:Filter> for my Attribute-Based SLD.

Here's what I've tried: I put my SLD in a text file and load it using SLD: some_url/text.sld param, instead of putting it straight to SLD_BODY like I did before. But another error happens, the map won't load (pink image) and I'm getting this error Resource interpreted as Image but transferred with MIME type application/vnd.ogc.se_xml: http://localhost/geoserver/wms?... Seems like the WMS request doesn't return an image but XML.


texts: How to map a point to a polygon label from a shapefile

texts:
I have a shapefile that breaks the city of San Francisco up into a number of different districts depending on who your supervisor is.

I need to make a lookup tool that takes a single long lat value and says which supervisor district the point lies in. The supervisor district is one of the attributes of the shapefile.

How would I go about doing this in an efficient manner?

Edit 1:
Seems like something like this is what I'm looking for. Trying to figure out how I access the attributes once an intersection is found


texts: Does the NED use geocentric or geodetic latitude?

texts:
What kind of latitude does the National Elevation Dataset (NED) use? Geocentric or geodetic or something else? I can't find the information anywhere on their site (http://ned.usgs.gov/index.html) or in the documentation that comes with the data I downloaded.

I'm new to GIS. If latitude type is unspecified, is there a standard default assumption?


texts:
I'm looking for historians' best estimate of the location of the Lewis & Clark expedition at each date during the trip, in a machine-readable format (e.g. that I could use in a GIS system).

(Since the party split up at various points, ideally I'd get the locations of each group on each date.)

After a bit of searching, I'm surprised not to have found anything like this. Can anyone point me in the right direction?


texts: What are some GIS outlets for my open source MapBasic utility module and library?

texts:
I'm actively developing a MapBasic module that provides additional MBX programming functionality by wrapping .NET 4.0 and Win32 API calls in easy to use MapBasic functions and subs.  MB_Utils acts as an abstraction layer so a user can easily do things that are not directly supported in MapBasic: get Windows user name, display custom progress bars and dialogs, perform file i/o, etc.  I've sort of been filling in the frustrating shortcomings of MapBasic I've come across.

I'm looking fro some suggestions on where to post MB_Utils so as to gain exposure to other MapBasic developers.  MB_Utils is hosted on GitHub:

https://github.com/mfcallahan/MB_Utils

https://github.com/mfcallahan/MB_Utils/releases


texts: Where do field aliases reside in an SDE database?

texts: If I want to get field aliases for an SDE table through pure SQL, what table do I query? My geodatabase was created with ArcGIS 10.1 in an Oracle database.


texts:
I follow this great post:pgRouting 2.0 for Windows quick guide, and I have loaded the test data from GeoJSON, like this:

But I don't know how to load this layer into my database, and my DB manager doesn't look the same as in Anita's post:

How should I do to load and test this layer into my PostGIS database?


texts: Sorting based on computed property (NDWI)

texts:
I am trying to sort a series of sentinel 2 images based on their average NDWI in the google earth engine. But nothing I do in sorting seems to work. Here I ll show what i tried:

's2_ndwi_mean' is a list of objects with 2 properties:

'NDWI' which is the average mean NDWI
'system:time_start' which is the time the image was taken.

This variable forms the starting point for this example.

Sorting a single dimensional array works:

for (var i = 0; i < 10; i++) {
    var ndwi_test = ee.Number(ee.Dictionary(s2_ndwi_mean.get(i)).get('NDWI')).float()
    ndwi_array.push(ndwi_test)
}

ndwi_array= ee.Array(ndwi_array).sort()

However sorting a list of objects based on their property does not seem to work. I tried a similar approach to this:

cars = [
    {
        name: "Honda",
        speed: 80
    },

    {
        name: "BMW",
        speed: 180
    },
]

cars.sort(function(a, b) {
    return a.speed - b.speed;
})

But I cannot get this to work for earth engine list types.

Then i tried sorting a two dimensional array instead of a list. This is how it should work for a Javascript Array. But this seemed impossible for an earth engine array:

var a = [[12, 'AAA'], [58, 'BBB'], [28, 'CCC'],[18, 'DDD']];

a.sort(sortFunction);

function sortFunction(a, b) {
    if (a[0] === b[0]) {
        return 0;
    }
    else {
        return (a[0] < b[0]) ? -1 : 1;
    }
}

Does anybody know what the approach should be. Sorting a list of images based on computed average NDWI's?


texts: Converting feature from draw interaction to WKT string

texts:
I've problems with converting the feature which I'm getting in draw interaction into the WKT string. I saw multiple examples which looks like it should be simple but I'm always getting the following exception while calling format.writeGeometry() function:

Here is my code:

function getDrawInteraction() {
    return new ol.interaction.Draw({
        source: this.source, //it is ol.source.Vector() for my layer
        stopClick: true,
        type: "Circle"
    });
}

function draw() {
    var drawInteraction = getDrawInteraction(),
        format = new ol.format.WKT();

    drawInteraction.on("drawend", function (event) {
        var feature = event.feature,
            wkt;

        wkt = format.writeGeometry(feature.getGeometry());

    });
    map.addInteraction(drawInteraction);
}

I'm able to draw a circle on the map but what I am doing wrong WKT string?


texts: Method to shade or overlay a raster map to reflect time of day and ambient light

texts:
I have a daytime raster of my study area.
I have 20+ counters in the study area that count "stuff" with a timestamp.
Using PyEphem I determined the type of ambient light for each collected data point.

What I'd like to do:

Add a grey/semi transparent overlay to reflect the ambient light (Dusk, Sunrise, day, sunset, and dark).
The end goal: produce a time-based visualization of the data collection that highlights the time of the day. Or be able to produce a snapshot at a certain time with the proper shaded overlay.
Think of this as this the image:

 source


texts: Retrieve last Row and specific field values with SearchCursor

texts:
I have a shapefile from which I want to retrieve the last row values for specific  fields to use them as input in a tool. So far this is what I have which does not work.

length=int(arcpy.GetCount_management(shapefile).getOutput(0))
print length
cursor=arcpy.SearchCursor(shapefile)
for row in cursor:
        row= row[length-1]
        endx= row.getValue("LastX")
        endy= row.getValue("LastY")
        sourcePoint= arcpy.Point(endx, endy)


texts:
After a few hours I am about to give up on finding a kind of documentation for the QGIS Python API that actually works/gives the correct information. It is obvious that there is no actual code documentation besides the one for the C++ API which differs from what works in Python. Ok, there is the help()-function, too, to get information about classes which is more than ust tedious to use. As an entry point, I thought I could use the tutorials/cookbook: http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/vector.html#modifying-vector-layers. Even simple things like adding features does not work the described way in the current QGIS Version 2.4.0-Chugiak (there is no addFeature-function for example).

I wonder how working with the Python console and writing plugins makes sense under such unfavorably conditions.

What sources can actually be used to lookup and learn the Python API for QGIS for current versions?


texts: Publish as a service from ArcGIS desktop to ArcGIS online without basemap

texts:
I'm trying to publish a map document as a service from ArcGIS Desktop 10.3.1 to ArcGIS Online.  It works just fine, except that ArcGIS Online automatically loads in a basemap.  This usually isn't an issue, but I'm working with a map of Io, one of Jupiter's moons.  A basemap of the earth is less than helpful.

I've seen the Javascript code to hide the basemap, but I was just curious if there's any way besides coding to remove the basemap.  I couldn't find a setting in the Publish a Service dialog box.

Is it possible to publish a map document as a service and only get the map document without ArcGIS Online's basemap?


texts: I have asked about this before.  If I import a USGS DRG raster for Anchorage, Alaska using EPSG:26706 - NAD27 / UTM zone 6N all the UTM coordinates are correct.  Now if I want to look at the coordinates in latitude/longitude and use ESPG:4267 NAD27 the image is rectified and stretched so the East/West looks twice as long as the North/South.  The latitude/longitude of the coordinate points are correct but the map just looks funny.  What am I doing wrong?


texts: Identify raster pixel value in ArcGIS Portal Map Viewer?

texts:
I have a raster I've published to ArcGIS Image Server 10.5.1. I'm trying to create a web map in 10.7.1 Portal where I can query the the pixel value. Ideally the user can click on the raster and a popup will display the pixel value.

Is this functionality possible? I tried publishing my raster as a Tiled Map Service and an Image Service with no luck.


texts: Cannot add plugins to QGIS on Ubuntu 11.10 Oneric

texts: I have installed Qgis on ubuntu. it works fine but I am not getting option to add new plugins only option to manage plugins is available. Please help


texts: Is there a way to determine the datum used when all you have is a lon and lat?

texts:
I have a dataset of US intersections from 2011 that I'm going to be using in a project together with data based in WGS84. Included in the CSV file is a variety of fields two of which are the lat and lon. The accompanying documentation makes no mention of the datum that was used for the coordinates. Is there a way to calculate it? If not, because it's US data, would it be safe to assume it's one based on NAD83?

I am going import this into a geography field in SQL Server 2012. I'm guessing that the SRID is either 4269 or 4326. I'd like to be able to determine which it is and transform it if necessary.


texts: ArcGIS 10.2: ArcPy - how to replace workspace paths?

texts:
I am teaching myself ArcPy and I am attempting to make a code to replace paths to the items in my workspace. I had a file geodatabase and its path changed so I would like to automatically modify the path with ArcPy.

>>> mxd=arcpy.mapping.MapDocument("CURRENT")
>>> arcpy.mapping.MapDocument.findAndReplaceWorkspacePaths(r"E:\GIS\_Wolverine Time\ArcGIS\California_DrinkingWater\Cali2.mdb",r"E:\GIS\California_DrinkingWater\Cali2.mdb")

I thought it was that simple but I received the following error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
TypeError: unbound method findAndReplaceWorkspacePaths() must be called with MapDocument instance as first argument (got str instance instead)


texts:
I'm trying to use the grass algorithm v.net.connect via the processing menu using a test set of 3 point objects and a single line object. All 4 objects only have geometry and an ID field.
With the parameters set to these two layers and the threshold still at 50.00 and the Network set to save to a temporary file, the algorithm appears to run but does not produce an output.
The log file appears ok bar a line near the end saying "ERROR: option : exists. to overwrite, use the --overwite flag"

My question is how do I set this flag to enable the output? or is this a bug with 2.18.13 as i have used this algorithm in the past on 2.14.


texts: Loop over columns of shapefile in modelbuilder

texts:
I use a shapefile (rain.shp) with a bunch of columns (rain_1979, rain_1990, rain_1991, ..., rain_2015). I would like to transform every column of the shapefile into a different raster using the Feature to Raster function in Modelbuilder (ArcGis 10.4.1 for Desktop). By manually selecting Field for one of the columns (e.g. rain_1989), everything works fine. But since I have so many columns, I'd prefer to use a loop.

I tried it with a 'For' loop from 1979 to 2015. My problem is however, that as soon as I try to use a placeholders in "Field" in the Feature to Raster function, it says that the Field rain_%Value% does not exist.

Is there any way to do this in Modelbuilder?


texts: ArcGIS Python Script Tool Field Name Drop Down?

texts:
I have a python script that checks file path names for corresponding files in a drive location. All of the inputs were originally hard coded, and it works well enough:

    import arcpy
    import os
    shp = r"C:\Users\dotwell\Desktop\Locations.shp"
    fields = ["ZdriveLoc"]
    f = open (r"C:\Users\dotwell\Desktop\ValMapsTest.txt", "w")
    cursor = arcpy.da.SearchCursor (shp, fields)
        for row in cursor:
            Roll=str(row[0])
            if not os.path.exists(Roll):
            f.write(Roll + os.linesep)

    f.close()

I'm trying to create an ArcGIS script tool from this, and I've only had partial success. I can use

shp = arcpy.GetParameterAsText(0)

to allow the user to select the shapefile. I'd like to then have a drop down list of the fields in that shapefile for the user to select the field that contains the path names, but I can't get it to work. Can someone show me how to generate the drop down? Here's my latest try. It works, but requires the user to type in the field name. It would also be nice for the user to designate the output text file name and location.

import arcpy
import os
shapefile = arcpy.GetParameterAsText(0)
field = arcpy.GetParameterAsText(1)
f = open (r"C:\Users\dotwell\Desktop\ValMapsTemp.txt", "w")
cursor = arcpy.da.SearchCursor (shapefile, field)
for row in cursor:
    Roll=str(row[0])
    if not os.path.exists(Roll):
        f.write(Roll + os.linesep)

f.close()

Here are the parameters in my ArcGIS script tool:


texts: Debugging "Hello world" sample plug-in : What about my breakpoint ...?

texts:
i'm trying to understand the Eclipse/Pydev/Qgis configuration :
i installed Eclipse (on a win7 pc) + Pydev and Remote Debug and the "hello world" example plugin in Qgis.

I open HelloWorld.py  in Eclipse and set a breakpoint at the beginning of the class definition. I start the debug server.

When i start Remote Debug in QGIS i get a "debuggin connection activated" message. all seems to be running ok ...

But when i run the Hello World plug-in in Qgis the execution is not stopped at breakpoint as expected and the message box "Hello World" pops up on the screen ...

What could i be missing ?


texts: Make PostgreSQL recover two geometry coordinates points split into columns

texts:
So I want to make a SQL that converts my PostgreSQL coordinates (x1,y1,x2,y2) to 4 different columns. By now this is what I have:

trim ('MULTILINESTRING((' from (split_part(ST_AsText(SSDMT.geom), ' ', 1))) as xi,
split_part(split_part(ST_AsText(SSDMT.geom), ' ', 2), ',', 1) as yi,
split_part(split_part(ST_AsText(SSDMT.geom), ' ', 2), ',', 2) as xf,
trim ('))' from (split_part(ST_AsText(SSDMT.geom), ' ', 3))) as yf

And it returns:

I also made this one:

trim('[[[' from (split_part((ST_AsGeoJSON(SSDMT.geom)::JSON->>'coordinates')::text, ',', 1))) as xi,
trim(']' from (split_part((ST_AsGeoJSON(SSDMT.geom)::JSON->>'coordinates')::text, ',', 2))) as yi,
trim('[' from (split_part((ST_AsGeoJSON(SSDMT.geom)::JSON->>'coordinates')::text, ',', 3))) as xf,
trim(']]]' from (split_part((ST_AsGeoJSON(SSDMT.geom)::JSON->>'coordinates')::text, ',', 4))) as yf

And it returns:

Is there a better way to do this?

The original query:

SELECT
SSDMT."COD_ID" ID_CABO,
SU."COD_ID" ID_SUB, SU."NOM" NOM_SUB,
CT."COD_ID" ID_ALIM, CT."NOM" NOM_ALIM,
SSDMT."FAS_CON" FASES,
SSDMT."COMP" COMPRIMENTO,
CBIT."COD_ID" COD_BITOLA,
CBIT."DESCR" NOME_BITOLA,
CMAT."COD_ID" COD_MATERIAL,
CMAT."DESCR" NOM_MATERIAL,
SC."R1",
SC."X1",
ST_AsGeoJSON(SSDMT.geom)::json->'coordinates' as COORD
from "BDGD2018"."SSDMT" SSDMT
join "BDGD2018"."SEGCON" SC on SSDMT."TIP_CND" = SC."COD_ID"
join "ANEEL"."TCABOBIT" CBIT on SC."BIT_FAS_1" = CBIT."COD_ID"
join "ANEEL"."TCABOMAT" CMAT on SC."MAT_FAS_1" = CMAT."COD_ID"
join "BDGD2018"."CTMT" CT on SSDMT."CTMT" = CT."COD_ID"
join "BDGD2018"."SUB" SU on SSDMT."SUB" = SU."COD_ID"
where SSDMT."CTMT" = 'NLAE733'

Returns:


texts: How to get all field names from table using arcpy.ListFields?

texts:
I have a table:

,

Script at ModelBuilder's "calculate values":

Expression:

a("%path%")

Code BLOCK:

def a(path):
  fields = arcpy.ListFields(path,"FID_*")
  for field in fields:
    return field.name

Gives me only ONE field name!


texts: How to delete a raster file after opening the raster in rioxarray

texts:
I wanted to able to delete a temporary file, after opening it with rioxarray and start to process the final product.
I followed two ideas:
a)
with rioxarray.open_rasterio(temp_file) as raster:
  # do processing
  # save output
os.remove(file)

b)
raster = rioxarray.open_rasterio(temp_file)
# do processing
# save result somewhere
del raster
gc.collect()
os.remove(file)

Both didn't work because, the file was still opened. How can this temporary file be deleted?


texts: osm2pgsql: does table planet_osm_roads contain any non-duplicate data?

texts:
Although there is a osm wiki page that attempts to shed light on the different tables which are created within a postgres/postgis database there is something that is not quite made clear:

Is the data in table planet_osm_roads duplicate? (i.e. would it provide anything that could not also be retrieved however indirectly from the other tables planet_osm_line,planet_osm_point,planet_osm_polygon ?

Bonus: If as assumed this is duplicate stuff is there a way to prevent its creation in the first place (assuming it increases the costs both in time and memory and disk space)?


texts:
I create a table with a geometry column (geom2)

 CREATE TABLE reu_oin_oi2 AS
 (SELECT eingangsdaten.*,
 ST_Intersection (test.geom, eingangsdaten.geom) as geom2
 FROM test, eingangsdaten
 JOIN test b
 ON ST_Intersects (eingangsdaten.geom, b.geom)
 WHERE test.id = 2);

Now some informations for the geometry column are missing (type of polygon ('Polygon') and SRID)
How can I add these inside the ST_Intersection command?


texts: PyQgis set Active Layer without using iface

texts:
I'm running my code in the linux console and there are always problems when using the iface function. I wanted to know if it is possible to define the active layer without using the iface function. Currently my code looks likes this:
qgis.utils.iface.setActiveLayer(rlayer)

rlayer = iface.activeLayer()

Is there another way to write this and get the same results?


texts:
While looking at this webpage, it seems there is an inconsistency for the time units handling.
Copying and pasting from the webpage:
from cftime import date2num
time_units = 'hours since {:%Y-%m-%d 00:00}'.format(times[0])
time_vals = date2num(times, time_units)

ds = xr.Dataset({'temperature': (['time', 'z', 'y', 'x'], temps, {'units':'Kelvin'})},
                 coords={'x_dist': (['x'], x, {'units':'km'}),
                         'y_dist': (['y'], y, {'units':'km'}),
                         'pressure': (['z'], press, {'units':'hPa'}),
                         'forecast_time': (['time'], times)
                })

ds.forecast_time.encoding['units'] = time_units

hence it seems that time_vals and not times should be assigned to forecast_time.


texts: How to make break values which are outside of the value range?

texts:
I want to compare mean temperatures for three future periods (2011-2040, 2041-2070, 2071-2100) with the measures mean temperatures of 1971-2000.
My resulting shapefiles will show the difference between the simulated and measured temperatures for each period.

I want to have a classification from -1°C. to +3°C. So I would like to have manual break values like -1,0, -0,8, -0,6 and so on. But how can I realize it? For example The difference between 1971-2000 and 2011-2040 is found within a range of -0,69 and +1,88. So the the "coldest" class begins at -0,69°C.

How can I create classes which are outside of the range of the values (for example -1,0 - -0,8)? When I try to edit the break values manually, it automatically jumps back to the coldest existing value (-0,69).


texts: Arcpy CopyRows writes only a fraction of MakeQueryLayer results

texts:
I have a query running against an Oracle non-spatial database that returns thousands of records, but when I run CopyRows against that query only 24 rows are carried over! What the heck is going on? I've tried with both ArcMap Python 2.7 and Pro Python 3.6.

query = '''SELECT * FROM ... WHERE ...'''
view = arcpy.management.MakeQueryLayer(
    input_database=input_database,
    out_layer_name=bph_layer,
    query=query,
    oid_fields="",
    shape_type="",
    srid="",
    spatial_reference="")
arcpy.management.CopyRows(view, bph_table)

print('Records from DB ({}): {}'.format(bph_layer, arcpy.management.GetCount(bph_layer)))
print('Records in View result ({}): {}'.format(view, arcpy.management.GetCount(view)))
print('Records out to {}: {}'.format(bph_table, arcpy.management.GetCount(bph_table)))

Emits:

Records from DB (BPH_plots_for_distribution_layer): 10139
Records in View result (BPH_plots_for_distribution_layer): 10139
Records out to BPH_plots_for_distribution: 24


texts: Suppressing parameter values from Python Window when running script in ArcGIS?

texts:
When I run a python script in ArcGIS Desktop 10.1 the values from the parameters will be written to the Python Window (marked with red in the picture below).

Does anybody know if it is possible to remove the parameter values or the other lines like the time stamp?


texts: How do I remove all clustering before add a new one

texts:
I would like to have a radio button to show the cluster which the first radio button: show clustering.The problem of my project is when I click at the first time, it shows clustering correctly, but when choose it for second time it shows again which clustering overlap each other (do not delete the first clustering)

    var command = L.control({position: 'topright'});
    command.onAdd = function(map){
    var div = L.DomUtil.create('div', 'command');
    div.innerHTML = '<input id="command" type="radio" name="abc">'
    document.getElementById ("command").addEventListener ("click",init,false)

    function init(){
          $.getJSON("XYpoint2.geojson",function(data){
          var clusters;
          var point = L.geoJson(data,{
              var marker1 = L.marker(latlng, {icon: myIcon1});
              return marker1;
          });
         clusters = L.markerClusterGroup();
            clusters.addLayer(point);
            map.addLayer(clusters);
    };

How can I delete the previous cluster when I click the radio button a second time?


texts: Are there Desktop GIS alternatives to ArcGIS 10.X for topology and vector conflation?

texts:
Is there any options other than ArcGIS for Desktop for topological rules and automatic topology integrity enforcing in a Desktop GIS? I am looking for something that will automatically snap features to some base layer bounds (eliminating gaps and overlaps in dataset).

I know about PostGIS topology functions but I would like something in a Desktop GIS.

QGIS 2 is trying to implement topological rules but nothing to automatically clean the data.

Is Esri's ArcGIS for Desktop actually the only solution?


texts: Dependence of point cloud data on camera parameters?

texts:
I want to do a 3d-reconstruction of a road surface using 'Structure from motion'. After taking the required image set of the road surface, I want to obtain the point cloud data using some open source applications (e.g. http://ccwu.me/vsfm/). Now , the question is regarding the spacing of the point cloud data that I will obtain after going through all the necessary steps i.e. :

How will the spacing between the data points depend upon my camera parameters? For example, is there a way to compare the spacing of the data points qualitatively or quantitatively if the image set is taken from say, a 15 MP camera and a 20 MP camera ?

i.e. I want to know comparatively how fine are my details of the road surface being captured from a , using the same example, if the image set is taken from a 15 MP camera and a 20 MP camera ?

Note : The data is of the form (x,y,z), where (x,y) refers to the location of the road surface and 'z' refers to the elevation of the road surface at that (x,y). Here , by spacing , I refer to the space between the location of the reconstructed road surface i.e. space between the (x,y) points .

Say, by using a 15 MP camera I got my data points, they have an average spacing of 2mm. Can I make a judgement beforehand that what will be my average spacing if I use a 20 MP camera and follow the same steps in the road surface reconstruction ? Can I estimate that it will be 1.5 mm (for example) ?


texts: How do I update joined data from MySQL programmatically?

texts:
I have a point vector layer (Sensores) in my QGIS project, I also have a MySQL database layer (SensoresDB) containing only data. I have joined these two layers in order to geographically reference my MySQL data and classify it.

My problem is the data in my vector layer (Sensores) doesn't update when the data in the MySQL layer (SensoresDB) does. That's ok, I don't need it to update in the exact moment it changes, but it would be nice if I could force it with code. I've tried using layer.setCacheImage(None) and layer.triggerRepaint() to no avail. Until now, the only way it does update is when I close the project and open it again.


texts: Simultaneous synchronisation of several QField files

texts:
I would like to synchronise several QField files using one command in QGIS. Is this possible with a Python script?
I would like to use the plugin QField sync but for the data of 10 tablets. So that I don't have to synchronise each tablet individually, a script would be very time-saving. Has anyone tried this before?


texts: PGRouting Layer plugin "UnboundLocalError: local variable 'db' referenced before assignment"

texts:
I am trying to do some basic computations using the PGRouting plugin for QGIS, but I am getting this error all time:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:/Users/Gonzalo/.qgis2/python/plugins\pgRoutingLayer\pgRoutingLayer.py", line 407, in run
    if db and db.con:
UnboundLocalError: local variable 'db' referenced before assignment

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.2.0-Valmiera Valmiera, c3a2817

Python path: ['C:/Users/Gonzalo/.qgis2/python/plugins\\processing', 'C:/PROGRA~1/QGISVA~1/apps/qgis/./python', u'C:/Users/Gonzalo/.qgis2/python', u'C:/Users/Gonzalo/.qgis2/python/plugins', 'C:/PROGRA~1/QGISVA~1/apps/qgis/./python/plugins', 'C:\\PROGRA~1\\QGISVA~1\\bin\\python27.zip', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\DLLs', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~1\\QGISVA~1\\bin', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\six-1.3.0-py2.7.egg', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\xlrd-0.9.2-py2.7.egg', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\xlwt-0.7.5-py2.7.egg', '.', 'C:\\PROGRA~1\\QGISVA~1\\apps\\qgis\\python\\plugins\\fTools\\tools']

I made a local installation of postgres 9.3.4, later I added postgis using the stack builder, and finally I followed the instructions posted here: PGRouting 2.0 for windows. I also installed python-psycopg2 for my version of python.

The SQL query runs well in pgAdminIII, but I can not make that QGIS plugin work.


texts: Not getting all decimal places when converting numbers to string (ArcObjects with VB .Net)

texts:
I would like to build some queries in my addin from layers in my map.  I am having problems with double numeric types. esriFieldTypeDouble has 14 - 18 decimal places or more depending on the database type. That's what I get when I build a "Select by Attributes" but when using ArcObjects with .Net I only get 12 or less decimal places (System.Collections.IEnumerator) and although I can get 14 decimal places when I retrieve a value directly from iRow.value() ,converting this number to a string in any way rounds it to 12.

I am testing this on SHAPE.STLength() and SHAPE.STArea() fields (SQLServer EXPRESS) and well as SHape_Area, SHape_Length (FGDB) and the problem is that since the value actually has more than those 12 decimal places I am getting in code but the actual value has more the query will not work.

Building the query in the attribute table's "Select by Attributes" correctly displays ALL decimal places in the Unique values list

Attempting to build the same list with the same ArcObjects that are uses to obtain this list "under the hood" System.Collections.IEnumerator rounds those to 12 places.

Is this a VB .Net issue?  Is there any way to get all the decimal places to string?


texts: Referencing external CSV file within ArcMap field calculator using python codeblock?

texts: I would like to know if there is a way to incorporate an external csv file into a custom made python script in ArcMap field calculator. I am inputting fields from the current layer, but need to access an external csv containing the mapping for these fields to convert them to values of interest. Is this at all possible?


texts: How to load a WMTS layer in QGIS using PyQGIS?

texts:
I've found tons of topics about how to load a WMS layer in PyQGIS, but nothing about WMTS.
In the GUI, adding a WMS or a WMTS layer is done by clicking the same button so I expect the syntax to be close, but I failed to figure it out entirely.

So I took this url (from this post) to make my testing :

http://opencache.statkart.no/gatekeeper/gk/gk.open_wmts?Version=1.0.0&service=wmts&request=getcapabilities

It works if you try to add layers from this source using the GUI.
I tried to copy the syntax that allows to add a WMS layer (so I added needed parameters like layers=, format=, style=) but my layer is never valid.

layer = QgsRasterLayer('url=http://opencache.statkart.no/gatekeeper/gk/gk.open_wmts?&layers=barentswatch_stedsnavn&Version=1.0.0&service=wmts&format=image/jpeg&styles=', 'anytitle', 'wms')
layer.isValid()
>> False

Does somebody see the mistake / know the trick (or where I can find documentation about it) ?


texts:
Using Postgresql 9.6, pgrouting 2.4.1  when I want to use pgr_dijkstra on a osm2pgrouting parsed database I get the following error:

ERROR:  Error, columns 'source', 'target' must be of type int4, 'cost' must be of type float8
CONTEXT:  PL/pgSQL function _pgr_parameter_check(text,text,boolean) line 40 at RAISE
PL/pgSQL function pgr_dijkstra(text,integer,integer,boolean,boolean) line 7 at assignment

I read that in previous versions this could happen, and tried to change the types alreay, using Qgis. I am not succeding.

What am I doing wrong?


texts: Difference in rounding of x and y point to pixel in raster

texts:
I am trying to round coordinates before querying a raster, to the precision of the raster, so that I can cache the results, and get the same values back even when the location is moved little. However, I find oddities in which pixel I access. In this case it's done through GDAL, but through Postgres it gives the same results.

gdalinfo ahn2_05m_dsm/r65dn1.tif:

Driver: GTiff/GeoTIFF
Files: ahn2_05m_dsm/r65dn1.tif
Size is 10000, 12500
Coordinate System is:
PROJCS["Amersfoort / RD New",
    GEOGCS["Amersfoort",
        DATUM["Amersfoort",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725],
            AUTHORITY["EPSG","6289"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4289"]],
    PROJECTION["Oblique_Stereographic"],
    PARAMETER["latitude_of_origin",52.15616055555555],
    PARAMETER["central_meridian",5.38763888888889],
    PARAMETER["scale_factor",0.9999079],
    PARAMETER["false_easting",155000],
    PARAMETER["false_northing",463000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","28992"]]
Origin = (30000.000000000000000,393750.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   30000.000,  393750.000) (  3d35'10.68"E, 51d31'11.12"N)
Lower Left  (   30000.000,  387500.000) (  3d35'18.69"E, 51d27'48.94"N)
Upper Right (   35000.000,  393750.000) (  3d39'29.96"E, 51d31'15.04"N)
Lower Right (   35000.000,  387500.000) (  3d39'37.64"E, 51d27'52.85"N)
Center      (   32500.000,  390625.000) (  3d37'24.25"E, 51d29'32.01"N)
Band 1 Block=10000x1 Type=Float32, ColorInterp=Gray
  NoData Value=-3.4028234663852886e+38

You can see below that when the boundary for the pixel in the X direction is on .99|.00, however, in the Y direction it is .00|01 (actually it is .000000000|.0000000001).

$ gdallocationinfo -geoloc ahn2_05m_dsm/r65dn1.tif 30714.999999999 391055.00
Report:
  Location: (1429P,5390L)
  Band 1:
    Value: 5.90999984741211
$ gdallocationinfo -geoloc ahn2_05m_dsm/r65dn1.tif 30715.00 391055.00
Report:
  Location: (1430P,5390L)
  Band 1:
    Value: 5.71999979019165
$ gdallocationinfo -geoloc ahn2_05m_dsm/r65dn1.tif 30715.00 391054.99
Report:
  Location: (1430P,5390L)
  Band 1:
    Value: 5.71999979019165
$ gdallocationinfo -geoloc ahn2_05m_dsm/r65dn1.tif 30715.00 391055.00
Report:
  Location: (1430P,5390L)
  Band 1:
    Value: 5.71999979019165
$ gdallocationinfo -geoloc ahn2_05m_dsm/r65dn1.tif 30715.00 391055.01
Report:
  Location: (1430P,5389L)
  Band 1:
    Value: 6.1100001335144

Why is this, and is there a way to correct it?


texts: Examining the covariance between biodiversity and other variables

texts:
I have som problems with an analysis and would be very grateful for any answers or suggestion.

The background is:
- I want to examine the covariance between ecosystem services in an area. I don´t really know if covariance is the correct word, however.

I have identified and combined different variables indicating high biodiversity using the weighted overlay and graded the areas from 1-10, whereas 10 is indicating high biodiversity.

Now I want to investigate how the biodiversity covariates/coexists/relates (english isn´t my native language) with other variables.

The other variables could be biotopes, soil type or other ecosystem services, like recreation.

My data consists of:
- polygons with the graded biodiversity and also the same areas in a raster.
- polygons indicating biotopes, soil type, precipitation etc.
- polygons indicating whether an area is used for a certain activity like mushrooming and so on.

The polygons are divided into green spaces at a county level.

So, what I basically want to do is to se if there are any hot spots where biodiversity are more likely to coexist with an other variable. I would like to visualize it on a map but also get a diagram over the relationships.

How do I do?

Oh, I noticed I wasnt allowed to post pictures. I hope you understand anyway!

Thanks in advance!


texts: Python Mapscript 6.0.1 layer.getShape API

texts:
I am attempting to follow the mapscript example described here - http://trac.osgeo.org/mapserver/wiki/RasterQuery

but ran into some trouble with

s = layer.getShape(result.shapeindex, result.tileindex)

It appears that this no longer works

In [63]: s = layer.getShape(result.shapeindex, result.tileindex)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/Users/calvin/work/luc/<ipython-input-63-969295303151> in <module>()
----> 1 s = layer.getShape(result.shapeindex, result.tileindex)

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/MapScript-6.0.1-py2.7-macosx-10.7-x86_64.egg/mapscript.pyc in getShape(self, *args)
   1395     def nextShape(self): return _mapscript.layerObj_nextShape(self)
   1396     def close(self): return _mapscript.layerObj_close(self)
-> 1397     def getShape(self, *args): return _mapscript.layerObj_getShape(self, *args)
   1398     def getNumResults(self): return _mapscript.layerObj_getNumResults(self)
   1399     def getResultsBounds(self): return _mapscript.layerObj_getResultsBounds(self)

TypeError: layerObj_getShape() takes exactly 2 arguments (3 given)

What is the argument I ought to supply to getShape method to get it working again?


texts: How fill the vias_path field for restriction table

texts:
In documentation explain the restriction table for pgr_trsp function, and that way I can restrict left turn or right turn {target_id, from_edge}. But if I want to restrict multiple links it doesn't explain how fill the via_path

In this case all roads are one way, and you can't do the red turns. So {2190, 2034, 673814} and {673814, 673816, 2190} are forbidden.

How I fill restriction table?

CREATE TABLE restrictions (
    rid BIGINT NOT NULL,
    to_cost FLOAT,
    target_id BIGINT,
    from_edge BIGINT,
    via_path TEXT
);

COPY restrictions (rid, to_cost, target_id, from_edge, via_path) FROM stdin WITH NULL '__NULL__' DELIMITER ',';
1,100,7,4,__NULL__
1,100,11,8,__NULL__
1,100,10,7,__NULL__
2,4,8,3,5
3,100,9,16,__NULL__
\.