texts: How to handle distance in postgis using ST_DWithin

texts:
I'm super new with this subject and was assigned the task of creating a query to find objects around a given point using some type of radius search using Postgresql with PostGis. After some research I came across the ST_DWithin method which is perfect for this.

I'm currently having some issues with the query. Here is how it looks like:

SELECT * FROM waypoints WHERE ST_DWithin(lonlat, '002000000100000f11c05e9335158b82804047dcae85b9e8c4', 10)

This is a table that stores lonlat using SRID=3857 and I'm creating the second geometry using  FACTORY.point(circle.fetch('lon'), circle.fetch('lat')) from Ruby RGEO

Now the documentation states that if I'm using geometries then the third parameter distance_of_srid should be in the units of the geometries.

The problem is that both lonlat and the point are using SRID 3857 which from what I understood correctly, uses units in meters. Logically I would expect that distance_of_srid would accept distances in meters. But when I send the query I get waypoints that are kilometers away not 10 meter away. But if I break it down to say .01 then it brings waypoints that are within 10 meters or at least that's what it seems.

I feel like I'm lacking some knowledge to overcome this but not sure what it is or where is a good place to get more acquainted with this subject. If someone could help me with this issue or point me to a better direction that would be awesome!

Update:

Based on a comment below I tried doing the following:

ST_DWithin(ST_Transform(archived_listings.lonlat, 4326)::geography, ST_GeogFromText('SRID=4326;POINT(-122.300115 47.724076)'), 10000.0)

But unfortunately I'm getting nothing with 10 km as before I was getting a house that was closer than that.


texts: Error clipping a raster with a shape layer "Attempt to create ADRG dataset with an illegal data type (Float32), only Byte supported by the format. "

texts:
I try to clip a raster with a vector layer. When I try to do it, I get the error "ERROR 1: Attempt to create ADRG dataset with an illegal data type (Float32), only Byte supported by the format."

Can someone explain what this means and how I could fix this?

Edit: additional information:

-im using qgis 2.18.9 (tagged).

-the file formats are .tiff and .shp (poly)

-I'm working on a win10 x64 system


texts: gdal_calc.py result are false and not the same as "direct" calculation

texts:
I need to do the raster calculator operation as "raster + constant" in Python. The input rasters are DEM data with float original type and -9999 for the NoData. The constant is float too.

In 1st way, I call the gdal_calc.py rutin:

command = 'gdal_calc.py -A {0} --outfile={1} --calc="A+{2}" \
       --NoDataValue=-9999 --type="Float32" --overwrite \
       --quiet'.format(filename, savename, median)
os.system(command)

In 2nd way, I open the file with gdal and directly operate the band values:

driver = gdal.GetDriverByName('GTIFF')
AD_file = gdal.Open(filename, gdal.GA_Update)
band=AD_file.GetRasterBand(1).ReadAsArray()
band = np.where([band > -9999], band + median, -9999) # raster calculator analog
AD_file.GetRasterBand(1).WriteArray(band[0,:,:])
AD_file=None

They both give the result. But... the results are different in therm of vales when I compare the min/max. For example, if the median=5 so the calculation is 'raster+5':

the source raster min/max: 241/642
the gdal_culc.py result min/max: 253/661
the 2nd way result min/max: 246/647

What happens with gdal_calc?!


texts: Converting time field omitting time and returning null values?

texts:
I have a time field in my shapefile that is formatted as HH:mm:ss.s
It is a text field, and I want to convert it to a date so I can calculate time elapsed across rows (if there's a better way to do this let me know).
Basically, each row is a point on a route, and I want to use the time field to calculate the amount of time elapsed between the first and last points on the route.

I tried using the convert time field tool, setting the input format as HH:mm:ss.s, and the output type as DATE. The output field however only shows values of 01/01/2001 for some rows and  for other rows.

time is the input field and time_c is the output field

I'm not sure what's going wrong here, here is what shows up in the geoprocessing results window:

I tried adding the date data to the field as well to make the format yyyy-MM-dd HH:mm:ss.s, but it still didn't work:


texts: How to automatically assign specific ColorBrewer ramp to polygons so that no touching polygons share the same color in QGIS 2.16

texts:
I'm creating a series of reference maps that need to differentiate between administrative regions within a country. I created a new color ramp using ColorBrewer Dark2 with 8 categories.

My administrative units are in sequential order, so when I classify them, the color ramp groups them in order, like so:

The polygons are also placed based on this order, so all of my administrative regions are roughly grouped like so:

This obviously doesn't work well--the point is to easily differentiate between different administrative units. I want each country to look more like this:

Thus far, I have been doing all stylizing by manually changing each administrative region's color so that none of the same colored polygons are touching and there are roughly the same number of each color within a country (so 3 pinks, 3 purples, 3 yellows, etc.). I want to be able to automate this process. I think that I need to use rule-based rendering, but I'm not sure how?

Can anyone point me in the right direction for automating this stylization? I get the feeling that I'll need to use some sort of if-then statement to loop through the regions, but this is not my area of expertise.

(Also, I'm using NaturalEarthData shapefiles with custom GIDs if anyone is interested in recreating my polygons.)


texts: Using shapely interpolate to evenly re-sample points on a linestring geodataframe at specified interval

texts:
I have a GeoPandas LineString GeoDataFrame. I am trying to evenly interpolate the spacing of the nodes that make up each line, with a couple constraints: keep endpoints unchanged, and specify maximum spacing interval. The original nodes should be removed if they are not needed (i.e., a node on a straight line versus a node at an angle on the line).
Desired output format is GeoPandas LineString GeoDataFrame with the following attributes: lineid, geometry.
I am using GeoPandas and Shapely libraries and would prefer to find a solution using Shapely. There is an arcpy tool that can do something similar (GeneratePointsAlongLines()), but I am looking for an open source solution.
I am trying to implement the solution here because it only requires Shapely (as opposed to solution here that requires OGR), but I am stuck on a couple parts:

my test script adds a column for every linestring. In the screenshot below, the first line should have nverts =22, and the second line should be 16 (in other words, the 0-column should read [22,16] and there should not be a 1-column)

when the return LineStr... part of my function is not commented out, I get an error: TypeError: 'Series' object cannot be interpreted as an integer I think this is maybe because it is not handling multiple lineids correctly.

Where am I going wrong in my script, or is there a better open-source, out-of-the-box solution for my line interpolation?
Using Windows 10, conda 4.8.2, Python 3.8.3
Libraries
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
from shapely import wkt
import pyproj
from pyproj import CRS

Setup repeatable example
##### BUILD EXAMPLE LINESTRING GEODATAFRAME #####
myid = [1, 1, 1, 2, 2, 2]
myorder = [1, 2, 3, 1, 2, 3]
x = [550338, 550428, 550523, 550529, 550689, 550703]
y = [3795929, 3795798, 3795659, 3795419, 3795290, 3795257]
myepsg = 32611

df = pd.DataFrame(
    list(zip(myid, myorder, x, y)),
    columns =['myid', 'myorder', 'x', 'y']
    )
gdf_pt = gpd.GeoDataFrame(
    df, geometry=gpd.points_from_xy(df['x'], df['y'])
    )
gdf_pt = gdf_pt.set_crs(epsg=myepsg)
gdf_line = (
    gdf_pt.sort_values(by='myorder')
    .groupby('myid')['geometry']
    .apply(lambda x: LineString(x.tolist()))
    )
gdf_line = gpd.GeoDataFrame(gdf_line, geometry='geometry')
gdf_line.reset_index(drop=False, inplace=True)

## delete stuff used to build because I don't
## have access to it in a real application
del myid, myorder, x, y, myepsg

Plotting
display(gdf_pt.style.hide_index())
display(gdf_line.style.hide_index())

ax = gdf_line.plot();
gdf_pt.plot(ax=ax)
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);
ax.xaxis.set_major_formatter(ScalarFormatter())
ax.ticklabel_format(style='plain', axis='both', useOffset=False)

## this is not necessary, but I like to add labels
for x,y,z1,z2 in zip(gdf_pt.x, gdf_pt.y, gdf_pt.myid, gdf_pt.myorder):
    label = str(int(z1)) + '-' + str(int(z2))
    plt.annotate(label, (x,y), textcoords = 'offset points',
                 xytext = (12,-5), ha = 'center')

Analysis & Problem Area
def mytest(row, gdf, maxlen):
    nverts = round(gdf.length / maxlen).astype(int)
    return nverts
    #return LineString(
        #[gdf.interpolate(float(n) / nverts, normalized=True)
        #for n in range(nverts + 1)])

lineid = gdf_line.myid
geometry = gdf_line.geometry

testout = (
    pd.DataFrame(list(zip(lineid, geometry)), columns =['lineid', 'geometry'])
    .pipe(gpd.GeoDataFrame, geometry=geometry)
    .pipe(lambda df: df.join(df.apply(mytest, gdf = gdf_line, maxlen = 15, axis=1)))
    )

display(testout)
print(type(testout))


texts: Producing singlepart featureclass with unique IDs for each spatially non-connected polygon?

texts:
What is the best method to convert a singlepart shapefile (i.e. one record that contains several non-connected polygons) into a mutlipart shapefile?

For clarification the original question asked how to produce a singlepart featureclass with unique IDs for each spatially non-connected polygons.

The explode method in Advanced Editing or Multipart to Singlepart would give this result.


texts: High Resolution Sinusoidal Interrupted Projection of the Earth

texts:
My grandfather built a rather large globe 10 years ago and never got to complete it (he started drawing about a fifth of it by hand and just gave up after that).

I have been looking for a sinusoidal interrupted projection, with a traditional linework (but any would do) that I could get printed and apply to his globe to finish his work for him (as a surprise!).

I have contacted a few companies but they are charging around 3000 EUR (linework and background topography image) which is way too much for me.

I am a bit lost with the science of it and would really appreciate if someone who understands that stuff could help me out of bit :
Basically where kind I find a high res image of the earth with country and city names, borders, etc... and how do I make it into a sinusoidal interrupted projection?

Data I got from my grandpa a couple years ago, when I was shopping around :
- globe circumference = 2 meters
- diameter = 0,63662m
- scale = 1/20 M
- 10° L=111km111111/20M
          =5cm555

Thanks!


texts: DN value from raster to vector tool (QGIS 2.0)

texts: I am using QGIS 2.0 to convert a DEM to vectors, then create polygons of zones of slopes. However, I am unsure of in what units my new "DN" column is. It has lost of precision of z-values from my original XYZ points, and am not sure how it compares.


texts: Openlayers 3 different projection Geojson view

texts:
My map projection is EPSG:3857. And I want to show a GeoJSON data which has EPSG:4326 projection.

            var source = new ol.source.Vector({
                url: 'http://localhost:15485/api/locations',
                format: new ol.format.GeoJSON()
            });
            var locationsGeojsonLayer = new ol.layer.Vector({
                source: source,
                projection: 'EPSG:3857'
            });

Locations are getting from url but does not apper on map. Openlayers does not automatically convert URL data.

if I use ajax call and add features manually it appears on map.

            $http({method: "GET", url: "http://localhost:15485/api/locations"}).then(function (response) {
                angular.forEach(response.data, function (item) {

                    var feature = (new ol.format.GeoJSON()).readFeature(item,{
                        dataProjection: 'EPSG:4326',
                        featureProjection: 'EPSG:3857'
                    });

                    source.addFeature(feature);
                })
            });


texts: ArcMap crashes on TabbedFeatureInspector sample

texts:
The TabbedFeatureInspector, (VB.NET), sample crashes ArcMap when I do the following:

Open ArcMap.
Open map document with feature class with class extension loaded.
Select a feature.
Start an edit session.
Click the Attributes button to open the Attribute editor.
Click the Create Features tab.

At this point the Create Features tab displays as a blank tab and ArcMap locks up.   The wierd thing is that this only happens with the initial selection of a feature and opening of the Attribute editor.  If you open the attribute editor, and then select a different feature before selecting the Create Feature tab ArcMap doesn't crash.

I initially found this by building a tabbed inspector into an application using the sample as the template.  After discovering the problem I tried it in the TabbedFeatureInspector sample and found the same result.

The offending code appears to the be SetParent line in Inspect event where the inspector parent window is set to the groupbox on the user control.

Public Sub Inspect(ByVal objects As ESRI.ArcGIS.Editor.IEnumRow, ByVal Editor As ESRI.ArcGIS.Editor.IEditor) Implements ESRI.ArcGIS.Editor.IObjectInspector.Inspect

    SetParent(Me.m_inspector.HWND, Me.GroupBox1.Handle.ToInt32)
    ShowWindow(Me.m_inspector.HWND, SW_SHOW)

Any suggestions on how I can correct this?


texts: How to zoom to feature when clicking on drop-down list box item with ArcGIS JS API?

texts:
I am developing a web application with ArcGIS JS Api that is similar to this smart campus: UCSB Interactive Campus Map. Here, in the Building drop down list box, you can select a building, and the service map next to the sidebar will simply zoom and select that polygon feature.

The thing is that in my application, I do have a sidebar with a drop down list box, where you can select buildings as well. But I don't know how to make the map service to zoom-and-select that building when the user clicks on a determined Building code, such as "Library" or "Engineering".

How to make this correspondence between a text item and a polygon?

Please, could you provide any ideas on how to make this web application behave the same as UCSB?


texts: Identify Coordinate Reference System in Tanzania?

texts:
I'm trying to georeference this map, but I don't recognise the CRS. It is for a region of Tanzania. Here is the image on imgur in case Stack Exchange compresses the image.

How would I go about IDing this myself in the future?


texts: Using Leaflet Js, is it possible to know the onclick location of a Marker (ignore anchor point)?

texts: I have a marker object that is of a specific height and width, and although it is pointing a single coordinate, the visual representation of the marker spans many pixels. When I click on the visual representation of the marker, I would like to get the underlying map coordinates, but instead I only have access to the single lat/lng coordinate associated with the marker.


texts: Popup marker info problems in Leaflet when using markercluster

texts:
I am new to Leaflet and recently did a course on using Leaflet to make maps. I followed the instructions to implement leaflet marker cluster, but, when trying to extract data to fill the popup from my GeoJSON, each popup contains the same data.
While the marker clustering works very well and picks up all the point sin the GeoJSON, the data going into there popup is always the same!
The code I have used looks like this:
var markers = L.markerClusterGroup({
  showCoverageOnHover: false,
});

json.features.forEach(function (feature) {
  markers
    .addLayer(
      L.marker([
        feature.geometry.coordinates[1],
        feature.geometry.coordinates[0],
      ])
    )
    .bindPopup(feature.properties.Name);
});
map.addLayer(markers);

How can I get the relevant info onto each popup?


texts: Does arcpy have a method for accessing related records within a cursor?

texts:
I have a line feature class that represent sampling segments.  I have a data table of field notes in the same Geodatabase that is connected to the feature class with a Relationship Class.

Is there a built in method for arcpy to access related records?  The work around would be to add another cursor for the second table, but you'd think that if you take the time to relate your data there would be a quick and easy method for accessing related records programatically.

Software: ArcGIS 10.1


texts: How to convert a geometry into LineString whose start and end point is equal?


texts: Postgis know if point is part of linestring

texts:
I have a linestring and I want to know if a point is in it, I have tried ST_Intersects but it only works if the point is exactly as in the linestring.

I want to find points which are in the linestring but are not defined in it.

Any thoughts on this?


texts: ReadAsArray leaves Qgs in unstable state

texts:
I have a QGIS plugin that reads a raster to get the elevation at vertex points in a polyline.  The issue is that even though a single call to ReadAsArray(...) appears to work, if I immediately exit QGIS, A Crash Dump dialog appears. What is wrong with the way I making the call?

def run(self):
# show the dialog
self.dlg.show()
# Run the dialog event loop
result = self.dlg.exec_()
# See if OK was pressed
if result:
    selectedVectorLayer = self.getSelectedVectorLayer()
    selectedRasterLayer = self.getSelectedRasterLayer()
    if self.validateSelectedLayers(selectedVectorLayer,selectedRasterLayer) is False:
        return;
    gdal.AllRegister()
    ds = gdal.Open(selectedRasterLayer.source(), GA_ReadOnly)
    if ds is None:
        QgsMessageLog.logMessage(selectedRasterLayer.source() + ' not found!', 'LineSlopes', QgsMessageLog.WARNING)
        return
    dataPoint = ds.GetRasterBand(1).ReadAsArray(122, 149, 1, 1)
    del ds


texts:
I'm trying to get the MODIS tile that I need, I can get it using online calculator http://landweb.nascom.nasa.gov/cgi-bin/developer/tilemap.cgi but I need to get it in code.

How I can convert lat, long to MOD09A1 H/V tile?


texts: Repeated GPS coordinates in data are creating stacked points on map in QGIS

texts:
I have 30+ events of data collection for each GPS location (previous person on project has simply added new rows for each data collection period). When I import the excel file (via CSV format), each location therefore has 30+ gps points stacked on top of one another.
How can I get a single GPS point to reference 30+ entries without also having 30+ repeated spots on the map?


texts: Hyperlink from document (Excel, Word, PDF) to a specified placemark inside the KML file

texts:
I have been searching on the web and I couldn't find a solution.

I would like to hyperlink from a cell in an Excel spreadsheet to a placemark in my .kml file.

If insert a link and attach a file straight like .../philipines.kml then I have all placemarks in Google Earth. I need one only. I know that I can extract every placemark to a separate .kml, but is there some quicker solution?

I tried to grip the placemark via id and then input into my Excel hyperlink like this:
.../philipines.kml#ninoy-aquino

Unfortunately the file opens the same like without this extension.


texts: Using Rasterio for calculating Volume in a raster file

texts:
I am calculating the volume in a raster .tif file.

My current approach is area * height where I calculate the area in the .tif file and multiply it by elevation value of .tif file.

I have used Raster Surface Volume algorithm of QGIS, Raster Volume of SAGA GIS and other techniques. But I want to use only Rasterio package.

Is there any in-built functionality of Rasterio which calculates the volume of a raster .tif file ?


texts: How to change the QGIS status message (bottom left) via qgis plugins

texts: I want to access and change the text which shows the current status of the process via python plugin.Can this be done and how?


texts: How I get the masking area data and export data in Google Earth Engine?

texts:
Here is my code in GEE:

var CO = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_CO')
         .filterDate('2018-07-04', '2018-07-05')
         .select("CO_column_number_density");

var CO_merged = CO.toList(CO.size());

var CO_1 = ee.Image(CO_merged.get(0));

var fireCount = ee.ImageCollection('MODIS/006/MOD14A1')
                .filterDate('2017-07-04').first();

var Fire = fireCount.select('FireMask').gte(8);

var CO_fire = CO_1.updateMask(Fire);

Map.setCenter(20,-8, 8);

Map.addLayer(CO_fire, {palette:'red'}, 'data1');

Here is my new layer:

My questions:

I want to get the red region (Masked by Fire) CO columns and the latitude and longitude information, finally export this red region data and geographic information about these red areas. so how can I get this information in GEE?


texts: How can I transform X Y coordinates of NAD83 Texas North Central to latitude and longitude, preferably in Python?

texts:
I am at my wit's end here!

For a project at work, I am plotting a simple map in Power BI that pinpoints some of our internal utility permitting activity. One of the datasets I have been given appears to have X and Y coordinates, and I was told that the coordinates are in NAD83 Texas North Central. The data set does not have latitude and longitude. I am looking for the python script to change them or even the formula.

Here are some examples of the coordinates (X, Y):

2389721.373 7120786.992
2404093.075 7121270.858
2399450.901 7108775.082
2390198.570 7122950.801
2388505.487 7129711.434

I have tried the following code, but I get lat and lon decimals that are nowhere near where they should be on the map.

conv = 0.3408

x1, y1 = 2389721.373 * conv, 7120786.992 * conv # feet to meters

p2 = Proj('4326', proj="utm", zone=14)
lon, lat = p2(x1, y1, inverse=True)

...which yields:

lat, lon
(21.91704812776455, -95.95674542629739)
# Nowhere near Denton, Texas

Any help would be immensely helpful. For reference, the lat lon points should end up in the Denton, Texas area.


texts: Is there a st_buffer inverse function, that returns a width estimation?

texts:
When you supply a geometry g1 and a width w  to the standard OGC function ST_Buffer,
it returns a new geometry g2,

  g2 = ST_Buffer(g1,w)

This question is about the inverse: a function where you can supply two related geometries — g1 and g2 satisfies ST_Within(g1,g2) — and returns the width parameter,

  w = Unbuffer(g1,g2)

Now the interpretation of w is "average width".
So, where is this function Unbuffer defined?  Are there references for it?  Is there a standard library for it?

Complementing after bounty.

Examples of applications (estimation of): street width, sidewalk width, river width, canal width, internal polygon semelhance, riparian vegetation buffer width, etc.

Example of proposals of credible sources:  Pavelsky & Smith 2008 ;  ArcMap Script Procedure.  They are not solutions here because are (very) CPU consumer, "river-specialized", and can not be implemented by OGC standards. They calculate average of river widthes by each of many perpendicular (of a central line) lines; so, they are "brute force" algorithms.

Example of proposal without credible and/or official sources: troubleshooter here. In this case, you  can get the bounty finding sources that makes credible the posted solution.


texts: Having ModelBuilder auto-detect work-space environment?

texts:
I have a model that I'm attempting to edit, that is used company wide. We have a structured folder template that all projects conform to. The model takes information from the personal geodatabase and returns a materials list in an excel file.

The tool is for non Arc users, and choosing file path for input and output locations is causing issues.

I need the model to automatically default to whatever database is open, and have it export the excel file to the corresponding folder location for that project.

I've been trying to figure this out on my own and had no luck.


texts:
I have a folder of 28 mxds that I'm trying to export all of them into 28 JPEGs. I'm having difficulty with the script and I'm hoping someone can help figure where I went wrong. Basically, I want to be export to run automatically, so that I don't have to open each mxd to export them into JPEGs.

    import arcpy, os
>>> folderPath = r"H:\Users\2015\Map15-0001\Zone Maps - for External Website"
>>> for filename in os.listdir(folderPath):
...     fullpath = os.path.join(folderPath, filename)
...     if os.path.isfile(fullpath):
...         basename, extension = os.path.splitext(fullpath)
...         if extension.lower() == ".mxd":
...             mxd = arcpy.mapping.MapDocument(fullpath)
...             df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
...             arcpy.mapping.ExportToJPEG(mxd, r"H:\Users\2015\Map15-0001\Zone Maps - for External Website\New"+basename+".jpg", df,resolution = 1000)

This is the error that I'm getting:

Runtime error
Traceback (most recent call last):
  File "<string>", line 8, in <module>
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\mapping.py", line 1026, in ExportToJPEG
    layout.exportToJPEG(*args)
AttributeError: DataFrameObject: Error in parsing arguments for ExportToJPEG

The mxds all have Data Driven Pages running on them. I was able to previously take all of the mxds and export separate PDFs.

I'm running ArcMap 10.3.


texts: How can I use ICommandBars interface to find my control on a custom toolbar in ArcMap

texts:
I have an ATL object (CMyCombo) which implements ICommand and IToolControl which contains a combobox on a toolbar.Its default interface is defined as

__interface IMyCombo
{
STDMETHOD(GetCurrentSelection(UINT* selectedIndex);
}

I want to find selected index of this combo box when I click on another command in same toolbar.I can find my ATL object using ICommandBars.Find method as a ICommandItem pointer,but when QI this pointer passing IMyCombo IID,it returns E_NOINTERFACE error,and if I use static_cast I can convert this pointer to CMyCombo* and I can see all local variables of this class but if I want to call GetCurrentSelection I receive another error.

may you help me to find how can I access this method inside my ATL object from another command?
Thanks in advance


texts: Errors when copy from shapefile to PostgreSQL using ogr2ogr

texts:
I am trying to upload data of MULTIPOINTS for listed buildings in Wales in shapefile format to a PostgreSQL database I have run the command below:
ogr2ogr -progress -update -append -f PostgreSQL "PG:dbname=dbtest2 host=localhost port=5432 user=user password=password" -lco LAUNDER=YES -lco GEOMETRY_NAME=geom -nln $LAYER -nlt MULTIPOINT -sql "select test_SHP_COLUMNS from test_SHP_LAYER" -dialect sqlite --config PG_USE_COPY YES -gt 655360 test_SHP

having in return this errors:
ERROR 1: COPY statement failed.
ERROR:  invalid byte sequence for encoding "UTF8": 0xe7 0x61 0x64
CONTEXT:  COPY listed_buildings, line 2

I have used the same code for a similar layer of Multipoints and it works perfectly, but some rows in the table of the layer I am trying to upload with the previous command.
Do you know how can I deal with this?


texts: Automate process of duplicating point features - ArcGIS 10.0

texts:
I need to learn a way in which I can normalize data and in the same process duplicate certain records of a point feature class. The issue is as follows. I am provided address data in the following format (a single field):

"1234, 1236, 1238 Travis Blvd"
or
"1234/ 1236/ 1238 Travis Blvd"

I currently manually select each individual record in this format, and then make a copy (Ctrl V, Control C) into the target layer (same as origin). In the case above, I would make two copies. I then manually edit the address to make each record unique. It is totally fine that the geography for all resulting points is identical. All other fields in the record should be preserved in this process.

Does anyone have a slick idea on how to automate (script) this task?


texts: Cannot convert csv file into point shapefile using ArcGIS for Desktop?

texts:
I can import my large csv file into ArcGis Desktop and I can view it in Catalog and everything looks fine. When I go to 'create feature class' it does not recognize the column that I need to use for the z value. It is not an option I can pick.

I'm using ArcGIS 10.2 for Desktop and trying to do this in Catalog.

This is what I get when I try to make the XY shapefile. No Z options.


texts: Map matching a stream of GPS points to possible bus routes & C# map matcher

texts:
I have streams of GPS points that are transmitted from a device aboard a bus in real time. I would like to use these points to infer which route the bus is following. I have routes stored as GPS point of bus stops, with a path shape which I have interpolated points every 20 meters. Currently I use the incoming points to narrow down the possible routes the vehicle is on using a finite state machine, by indicating a vehicle as possibly being on a route if it is within 50m of any point on the route.

The value of 50m is required due to the inaccuracy of the GPS unit. I therefore figured it would be better if I used map-matching to snap the incoming GPS coordinate to a road first so that I could use a smaller value of say 10m. It then occurred to me that instead of matching to a general road network, I could match to a map of just bus routes.

In the case of bus routes however, several routes follow the same road. I figured I could thus possibly use a map-matcher to solve my whole problem.

Can any map-matching algorithm match to more than one road simultaneously?
I am writing my algorithm in C#. Does anyone know of any open-source C# implementation of a map-matcher I can tinker with?


texts: Per row commit when editing PostGIS data with QGIS, with PostGIS triggers checking inserted data

texts:
I have a PostGIS trigger assuring that records with a Point-geometry inserted in a table fall within Polygons in another table.

Simple example:

-- Trigger function
CREATE OR REPLACE FUNCTION check_pointWithinPolygon()
RETURNS trigger AS $body$
BEGIN
  IF TG_OP = 'INSERT'
  THEN
    IF ((SELECT count(*) FROM (SELECT poly.id FROM poly WHERE st_Within(NEW.geom, poly.geom)) as foo) > 0) THEN
        RETURN NEW;
    ELSE
        RAISE EXCEPTION 'inserted feature ID:% is not within a poly', NEW.id;
    END IF;
  END IF;
END;
$body$ LANGUAGE plpgsql;
COMMIT;

-- apply the trigger function in a trigger
CREATE TRIGGER myPointInPoly BEFORE INSERT OR UPDATE
ON points FOR EACH ROW
EXECUTE PROCEDURE check_pointWithinPolygon();
COMMIT;

All this works well. I am using QGIS to edit the layer. Some point fall within polygons ( ok) some outside ( and are caught by the trigger). The problem comes when QGIS "save edits" committs the changes.

Currently, all edits are discarded if even a single record is not accepted. I would like to commit per row, so that those that are correct, are inserted. I am not sure whether this is something that can be caught in the trigger logic ( likely), or also needs reflection in the QGIS commit code.


texts:
I need to find the longest string of ALL cells in a table.
Unfortunately, the table can have varying number of fields each time.
This need to be a generic function. Any leads please?


texts: How to open a geodatabase stored on a local server using ArcGIS runtime version 100?

texts:
I want to open a geodatabase stored on a local server.  I need to open this geodatabase so I can use its feature tables but I have not been able to find how to do it?

I am using c# .net and arcgisruntime v 100 , I have 3D cone graphics displayed on my scene, and I have an SQL server geodatabase called "Test".

Connection is established via arcCatalog, so I have an .sde file that stores my connection to my geodatabase localised on our local server.

I wanted to open this geodatabase in order to load its feature table so I can add features on to it. I used the following code:

code private async void addFeature(MapPoint structureLocation , string type , string cause) { var gdb = await Esri.ArcGISRuntime.Data.Geodatabase.OpenAsync(@"C:\Users\Manel\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\Connection to frparsynoserver.sde"); var table = gdb.GeodatabaseFeatureTables.FirstOrDefault(t => t.TableName == "Test"); await table.LoadAsync(); if (table.LoadStatus == Esri.ArcGISRuntime.LoadStatus.Loaded) {
        //  var lyr = new FeatureLayer(table);
        //  MySceneView.Scene.OperationalLayers.Add(lyr);

        var attributes = new Dictionary<string, object>();
        attributes.Add("typdamage", type);
        attributes.Add("primcause", cause);

        //create a new feature
        var newFeature = table.CreateFeature(attributes, structureLocation);
        await table.AddFeatureAsync(newFeature);

    }

}

But it didn't work?

How to open a geodtabase using .sde file connection with arcgisruntime version 100?

Can any one tell me how to do that?


texts:
Using ArcGIS 10.3 Desktop.

We have a single feature class within our geodatabase which is a merge of all the other feature classes in the geodatabase.  I would like to build a model that we can run everytime we input new data which will automatically take ONLY the new data that hasn't yet been merged in the other feature classes (and ideally look for and replace any of the older data that has had an attribute changed) and merge those values with the existing "AllFeatures" feature class.

Is this possible?  What tool(s) would I use?


texts: Why does imported CSV string data have length and precision of 0, and how do I change it?

texts:
When I import a CSV file of census data, the length and precision of all of the fields are zero.

As a result, when I try to join the census data with a TIGER file of corresponding census block groups, the joined fields are "NULL."

How do I prepare the CSV file for an effective join? Both the data of the CSV file and the field on which the join is based in the TIGER file are string format. Thank you for your help!


texts: Zooming to selected features with ArcObjects

texts:
I am trying to figure out how to set the active view extent to be able to view all the selected features on the ArcMap.

The Map has only one layer.
The features are filtered with attributes.

Below is the code where I tried,stuck with the part where I can loop through the features in
the feature selection set but not able to set the active view extent to zoom to show all the selected ones.

This is similar to doing the right click on the attributes table and selecting multiple rows and doing a Zoom to selected.

private void ZoomToMultipleDGVSelection(List<int> selectedRightIDs)
    {

        IMxDocument pMxDoc = ArcMap.Document;
        IMap pMap = (IMap)pMxDoc.ActiveView;

        ESRI.ArcGIS.Carto.ILayer layer = GetLayersClass.GetFieldBoundaryLayer;
        if (layer is ESRI.ArcGIS.Carto.IGroupLayer)
        {

            ESRI.ArcGIS.Carto.IGroupLayer groupLayer = layer as    ESRI.ArcGIS.Carto.IGroupLayer;
            ICompositeLayer pCompositeLayer = layer as ICompositeLayer;
            int layers = pCompositeLayer.Count;
            ILayer pLayer = pCompositeLayer.Layer[0];
            IFeatureLayer pFeatureLayer = (IFeatureLayer)pLayer;

            IFeatureSelection pFeatureSelection = (IFeatureSelection)pLayer;
            ISelectionSet pSelectionSet = pFeatureSelection.SelectionSet;

            IFeatureClass pFeatureClass = pFeatureLayer.FeatureClass;
            IQueryFilter pFilter = new QueryFilterClass();

            foreach(int ID  in selectedRightIDs)
            {
            pFilter.WhereClause = "RightID = " + ID.ToString();

            IFeatureCursor pFeatureCursor = pFeatureClass.Search(pFilter, false);
            IFeature pFeature = pFeatureCursor.NextFeature();

            pFeatureSelection.Add(pFeature);

            }

            //if (pFeature == null)
            //{
            //    System.Windows.Forms.MessageBox.Show("This section doesn't exist");
            //    return;
            //}

            IGeometry pgeom = (IGeometry)pFeature.Shape;
            pMap.SelectByShape(pgeom, null, false);
            IEnvelope pEnv = pgeom.Envelope;
            pMxDoc.ActiveView.Extent = pEnv;
            pMxDoc.ActiveView.Refresh();
        }
    }


texts: How to create a surface layer from gridded point data of air temperature and Precipitation?

texts: I want to aggregate monthly averages of temperature and precipitation over a 10-year period into weighted seasonal averages (DJF) for use in interpolation. I have downloaded Shapefiles for each month (i.e. DJF) and each year from NCAR (http://gisclimatechange.ucar.edu/). One way to do this would be to copy-paste the attribute table values and average that way but I'm sure there is a more efficient method. I'm wondering if QGIS has any sort of plugin or preexisting tool which could make the process easier.


texts: Create overlapping polygons from contour lines

texts:
Working in ArcGIS 10, I have a set of contour lines representing density of commercial square footage, and I have a set of points representing local maxima for each peak of the dataset. For each maximum I would like to find the largest contour line that contains only that maximum and no others.

My strategy is to convert the contours to polygons, spatially join the polygons to the points, choose only the polygons that overlapped 1 point, and for each maximum choose the largest polygon, an operation that will probably happen in MS Access.

I've run into one problem - when I convert the contours to polygons using "feature to polygon", the polygons don't overlap, so my maxima only join to one polygon, breaking the method.

How can I convert polylines to overlapping polygons? If possible, I'd prefer a solution I can fold into an ArcPy script.


texts: QGIS3 call model from windows path with Python script

texts:
I was successful in triggering a model with the Python editor with this one-liner:

    processing.runAndLoadResults('model:my_model', {'param1':sParam1,'param2': sParam2,'coordinatesystem':QgsCoordinateReferenceSystem('EPSG:25832')),'qgis:createpointslayerfromtable_1:Trees':'memory:'})

It only works, when the model is already added to the QGIS project. In my scenario this is not the case, which is why I need to replace the part 'model:my_model' with a windows path to the model.

I found some helpful info here:
Using processing.runandload in QGIS 3
The example is unfortunately for QgsRasterLayer . I assumed something like
"QgsModel" would serve my purpose, so I searched for a class containing "model" in the QGIS API Documentation, but without success.

How do you use processing.runAndLoadResults() with a model stored on disk?


texts: Raster Analysis - Identifying features from Raster

texts:
I have an IDW, tension spline and hill-shade rasters that show an archaeological landscape and from these rasters I want to identify new features which are seen in the landscape as depressions. And from this raster I want to calculate the size and shape of these features. I have not yet come across a form of analysis that would do this, or I am just not understanding the different types of analysis I have looked at.

I have been using GRASS, QGis, ArcMap/Secne and surfer, but I'm willing to give anything a try.

Thanks

Christina


texts:
I am trying to create a single side distance buffer around a line on QGIS.

Here is the code I am writing but I don't get the expected result:

processing.runalg("qgis:fixeddistancebuffer","input.shp",0.5821,25,False,"output.shp",side=left)

Any idea?


texts:
Similar to the provided example of selecting different layers in CartoVL, is it possible to select different basemaps?

For example, my map uses:

const map = new mapboxgl.Map({
  container: 'map',
  style: carto.basemaps.voyager,
  center: [24.5,-28.6],
  zoom: 4.85,
  scrollZoom: false,
});

And I would like to change from the voyager to darkmatter basemap when a function is triggered using a radio button.


texts: Converted Mapinfo file (originally in Chinese) to SHP but the attributes are jumbled/wrong

texts:
I converted a set of MapInfo files to SHP using org2gui. Unfortunately, the attributes got jumbled. The original dataset was in Chinese, so that is probably contributing to the jumbling. Screenshot below shows jumbles.

I uploaded the source Mapinfo files to filebin here, and I uploaded the converted SHP to filebin here.

Is there anyway to recover the jumbled attributes? I have opened SHP in ArcGIS with Chinese characters in the attributes before, so it's not solely a language character rendering issue.


texts: Setting label weighting for legend in ArcMap?

texts:
I am using data driven pages to produce a series of maps. For each map extent, surrounding areas are labelled. My legend is located in the corner of the data frame (realizing now that that wasn't a good idea).

Is there a way to ensure that the labels in my map won't be covered by the legend box, like assigning a weight to it?


texts: arcpy.CheckExtension() Works in ArcMap but not in PyScripter

texts:
I've had no trouble with using similar licensing in Python scripts previously, so I'm wondering if a recent upgrade/replacement of my PyScripter environment may have caused it to stop working.

I'm using ArcGIS 10.7.1.

Running the following in ArcMap works:

>>> arcpy.CheckExtension("Spatial")
u'Available'

But running the following in PyScripter (3.6.1):

import arcpy
print arcpy.CheckExtension("Spatial")

it just prints:

NotLicensed

How can I get this to work in the upgraded PyScripter like it used to in the old PyScripter environment?

NB:  If I run it on a PC where PyScripter has NOT been upgraded/replaced, it does work, and prints:

Available


texts: How to open SMAP's Surface Moisture data files

texts:
I downloaded surface moisture data from Nasa's Soil Moisture Active Passive (SMAP) sensor website (https://nsidc.org/data/smap/smap-data.html).

For each processing level (1-4), three types of files are present inside each folder: a *.qa file, an *.xml file, and a *.h5file.

I tried to open the h5 file with ENVI but the format could not be recognized on the Dataset Browser.

Gdal also could not read it as:

dataset = gdal.Open('SMAP_L1C_S0_HIRES_02045_D_20150620T094355_R11850_001.h5')
print dataset.RasterCount

returned 0 which means that the bands could not be read.


texts: Creating a mosaic out of rotated raster with gdal_translate

texts:
I have a number of images (the amount does not really matter) which are georeferenced with a correct transform ( since they show correct position and rotation in qgis).
But when I try to use gdal_translate the function struggles to output a correct mosaic if the images are rotated more than 2 degrees from the North.
It seems like the function has problems creating correct sized windows for each rotated image. I am wondering if anyone has stumbled upon this problem before and has a solution?


texts: Unable to debug Esri Addin for ArcCatalog

texts:
It was a toss of a coin whether this question is a GIS.SE or StackOverflow question... I'm going with GIS.SE because the problem seems to lie with ArcCatalog.

I have several addins that have been written for ArcCatalog (and also ArcMap); previously I was able to F5 debug the code using breakpoints now all my breakpoints are inactive with the message The breakpoint will not currently be hit. No symbols have been loaded for this document:

I have tried the suggestions from Stack Overflow:

Fixing “The breakpoint will not currently be hit. No symbols have been loaded for this document.”

This Breakpoint Will Not Currently Be Hit, The Source Code is Different from the Original (even after cleaning/rebuilding)

The breakpoint will not currently be hit. No symbols have been loaded for this document

and also

Why can't the breakpoint be "hit" when debugging an ArcGIS 10 Add-In? but this appears to be a different problem.

As this is an Addin (dll) so there is no startup project. I have cleansed all other versions. I made a subtle change to the Esriaddinx file in the message= to confirm the changes were being made - all good, there is no older version.

Whilst cleansing I noticed that addins for ArcCatalog are copied to %USERPROFILE%\Documents\ArcGIS\AddIns\Desktop10.X\{AddInID}\ where ArcMap addins are not.. I have tried copying the pdb file to the addin folder but it seems that the code is no longer being debugged and ArcCatalog has switched/redirected to the copy.

Is there a way to debug with breakpoints addins for ArcCatalog? Note that ArcMap addins and standalone applications debug just fine.

Environment:

Windows 7 64bit
ArcGis 10.1
Visual Studio 2010 Professional (both VB.net and C#)


texts: Get wind direction from raster stack U and V in R

texts:
I have this NetCDF file (about 130 kb size). I want to get for each day (layer) the wind direction (in meteorological or geographical directions) at every cell.

https://drive.google.com/open?id=1UxYGJK1-gIgXWZT6tE10bavmiNEpsVx3

I tried:

library(raster)

u <- brick("uv_comp.nc", varname = "u10")
v <- brick("uv_comp.nc", varname = "v10")

dirs <- atan2(v[[1]], u[[1]]) / (2*3.14159265) + (0.5 * 360) # for the first layer

However, I am not very sure that this method is the correct one. Can anyone confirm, please?


texts: Clipping Modis ocean color chlor-a data using shapefile in R?

texts:
How do I clip Modis ocean color chlor-a data using a shapefile in R?

.libPaths("C:/R-packages")
library(raster)
library(ncdf4)
setwd("G:/Geoinformatics/Second semester/Spatial Statistics/Spatial_Statistics/Practicals/Practical_3/MOD28_SST_daytime_8day_2017")
getwd()
dir()
#open .nc
nc_open ("G:/Geoinformatics/Second semester/Spatial Statistics/Spatial_Statistics/Practicals/Practical_3/MOD28_SST_daytime_8day_2017/A20170012017008.L3m_8D_SST_sst_4km.nc")

#projection

r2 <- raster("G:/Geoinformatics/Second semester/Spatial Statistics/Spatial_Statistics/Practicals/Practical_3/MOD28_SST_daytime_8day_2017/A20170012017008.L3m_8D_SST_sst_4km.nc",var="sst")
proj4string(r2)


texts: Getting QGIS installation path in PyQGIS?

texts:
Is there any way to get the QGIS installation path programmatically?

I'm working on something similar to plugin installer and it would be great to have it.


texts: Using the mouse to select records in ArcPy

texts:
I am running ArcGIS 10.1 Basic license with Python 2.7.2.

I have a two spatial files: the first listing vessel positions and the second activity on board.  I also have a collection of reports in Excel which list events I must investigate by looking at that spatial data.  All three files include an ID and a date-time field.

I want to be able to view the spreadsheet, using the mouse to work my way through the spreadsheet selecting events and then using the ID and the date-time to view the corresponding spatial data.  I am thinking a Definition Query and then Zoom to Layer.

I need the flexibility to jump around in the spreadsheet -- perhaps weeks later jumping back into the spreadsheet and following up on the original investigation of a particular event.  So, cycling through the spreadsheet in order with a cursor is not what I want.

I would also like to be able to update a field listing the status of the investigation.

The challenge, of course, is accessing that record manually selected from the spreadsheet and pulling out the ID and date time and updating the status of the event.

I have considered a few of approaches.

Select the record outside of Python and then use button to activate a tool.  The problem is that one cannot access the selection with Python.
Start in the Python script and then select the record.  Does anyone know how to display a table view in Python and then make a selection manually?
Select the record manually outside of Python and then export the record to a named table or file.  Start a Python script that always works with the named table or file and its single record.

I have searched this site and have not found anyone trying to do this.  Can anyone offer advice?


texts: Can't figure out WFS URL filter - what amI doing wrong?

texts:
I'm trying to figure out how to use the filtering function on a WFS layer. Any where I look, there is not an example of the XML AND the query alongside it explaining how it works. In the URL query below, my filter is completely ignored. (This is just an example, I have tried several modifications to this URL, to no avail.)

Can someone provide an example of a working URL query with the filter parameter?

http://maps.co.pueblo.co.us/outside/services/pueblo_county_roads_wms_only/MapServer/WFSServer?request=GetFeature&service=WFS&typename=pueblo_county_roads_wms_only:Roads&maxfeatures=2&filter=%3CPropertyIsLike%3E%3CPropertyName%3EOBJECTID%3C/PropertyName%3E%3CLiteral%3E66580%3C/Literal%3E%3C/PropertyIsLike%3E&


texts: I am new in working with GeoServer. At the beginning I tried to publish a WFS service. This works great, but this service consists of huge data. That's why I want it to start displaying the polygons at a scale of 1:2500. Is there a way to set this option in GeoServer?


texts: Creating a reference line among LiDAR points and retrieving distances between them?

texts:
I have a LAS dataset which contains points with x and y coordinates. For further analysis I need a new dimension which can be adjusted as the dataset requires. To make it understandable I made this (beautiful paint drawing):

So in this case I need a V values which shows the distance along the V (North-East; South-West)

How is it possible to create these values?

I'd like to do it in Python, but any other working method is welcomed.

(I have ArcGIS 10.2.2 [built in 2.7 python], python 3.5)


texts: Plot times series for stack raster based on a polygons in my area

texts:
I am using R language and I would like to create a EVI MODIS time series plot, EVI vs time. I can did  method to create a raster stack and plot the average EVI of the total area over time.  But I would like to plot for a select polygons : 140 polygons with centroid or lat/long coordinate for centroid  for polygon) within the area.

i have 1 raster for one month = 8 raster  for year :

These are the polygons link for polygon

https://drive.google.com/drive/u/0/folders/1_frvUnKXF1ASD4c65pHs-asUXJUveyh2

This is my script for plot avrage EVI=meanEVI of total area over times ( EVI/times) :

 rm(list=ls())
library(scales)
library(raster)
library(rgdal)
library(ggplot2)
library(rts)
library(sf)

# Create list of NDVI file paths
all_EVI <- list.files("D:/Rteledetection/Pivots/Doy2000",full.names = TRUE, pattern = ".tif$")
# Create a time series raster stack
EVI_stack <- stack(all_EVI)
# apply scale factor
EVI_EVI_stack <- EVI_stack/10000
EVI_EVI_stack
plot(EVI_EVI_stack)
summary(EVI_EVI_stack)

#Calculate Average NDVI
# calculate mean NDVI for each raster
avg_EVI_stack <- cellStats(EVI_EVI_stack,mean)
# convert output array to data.frame
avg_EVI_stack <- as.data.frame(avg_EVI_stack)

# view column name slot
names(avg_EVI_stack)

#renomer le nom avec meanEVI
names(avg_EVI_stack) <- "meanEVI"

# view cleaned column names
names(avg_EVI_stack)

avg_EVI_stack$product <- "MOD13Q1"

# add a "year" column to our data
avg_EVI_stack$year <- "2000"

julianDays <- gsub(pattern = "MOD13Q1.006__250m_16_days_EVI_doy2000|_aid0001","", row.names(avg_EVI_stack))

head(julianDays)

# add julianDay values as a column in the data frame
avg_EVI_stack$julianDay <- julianDays
# set the origin for the julian date (1 Jan 2000)
origin <- as.Date("2000-01-01")
# convert "julianDay" from class character to integer
avg_EVI_stack$julianDay <- as.integer(avg_EVI_stack$julianDay)
avg_EVI_stack$Date<- origin + (avg_EVI_stack$julianDay-1)

#if run!
head(avg_EVI_stack$Date)
# What are the classes of the two columns now?
class(avg_EVI_stack$Date)

# plot EVI
f <- ggplot(avg_EVI_stack, aes(x= Date, y=meanEVI), na.rm=TRUE) +
  geom_line(color = "blue", size = 1)+
  geom_point(size=2,color="black",pch=1)+
  ggtitle("MODIS EVI - 2000 Mato Grosso  Site") +
  xlab("Date") + ylab("Mean EVI_No units (Scale Factor = 0.0001)")+
  theme(text = element_text(size=10))+theme_classic()+
  scale_x_date(labels = date_format("%m-%Y"))+
  scale_y_continuous(name="Mean EVI_No units (Scale Factor = 0.0001)", limits=c(0.3, 0.5))+
  theme_bw()

plot(f)

and these results:

I'm trying to calculate the centroid for my polygons like this :

 #calcul centroid

    centroid = gCentroid(polyk,byid=TRUE)
    plot(polyk)
    points(coordinates(polyk),pch=1)
    points(centroid,pch=2)

but is not a good idea i need a simple method to plot times series for my polygon within the area. I know I can us extract::raster but I don't know how and where can I write this line of script !?
i have idea if i do vector to raster to my polygons and i do mask and i multiply the mask with my raster stack and calculate sum for all raster stack after and plot it but i don't know how i can do this script. Can you help with an explanation and example of what I need?

I know my script is not good but I am trying.


texts: Does Y mean latitude and X mean longitude in every GIS software?

texts: I am using Mapinfo and it has Y as latitude and X as longitude. Is that the same case for all mapping software? As for any country their respective value is multiple of 1 or -1. So for Nepal can I say it is on positive side +1 for both latitude and longitude? And for USA to be +1 Y and -1 X.


texts: Joining attributes of two line feature classes (street network, ArcGIS)

texts:
I have two line feature classes for a street network.  Their geometries are only crudely similar, while their attributes are different.  I am trying to add the attributes of one feature class to the other based on their closest line segments.  The screenshot is an example.  I want to assign the attributes of the middle red line to the selected lines of the the other feature class (and ultimately do it for the whole feature class).  Its been accomplished before, but I don't know how.

Any ideas or suggestions?


texts:
This is my code to draw a polygon on map but instead on putting a circle on each point, I want to put a image (custom icon). How do I do that?

var draw = new MapboxDraw({//modifier des buttons pour dessiner
    displayControlsDefault: false,
    controls: {
        line_string: true,
        polygon:true,
        trash: true,
    },
    userProperties: true,
      styles: [
          {
          "id": "gl-draw-polygon-stroke-active",//change le style de ligne de polygon
          "type": "line",

          "layout": {
            "line-cap": "round",
            "line-join": "round"
          },
          "paint": {
            "line-color": "#D20C0C",
            "line-dasharray": [0.2, 2],
            "line-width": 2
          }
        },
        {
          "id": "gl-draw-polygon-and-line-vertex-active",//change le style des points dd polygon
          "type": "circle",
          "filter": ["all", ["==", "meta", "vertex"], ["==", "$type", "Point"], ["!=", "mode", "static"]],
          "paint": {
            "circle-radius": 5,
            "circle-color": "#3c33ff",
          }
        }
       ]
    });


texts:
I'm using ArcGIS 9.3.

I've downloaded RFE (Rainfall Estimate) images from FEWS NET, and I'm trying to load them on ArcMap, to clip them using Mozambique's shape. But as the RFE images related to Southern Africa come with unknown coordinate system, I'm having problems projecting RFE and Mozambique shape on the same map.

I also tried to put them both at the same GCS which is WGS_1984, and then use the same projection for both. It didn't work.

I'm a GIS beginner. Which is my mistake in this projecting and reprojecting issue?


texts:
I'm using QGIS 2.18.9 on Windows 10. I've created an atlas in composer and when exporting images, the first few will export fine and then the map image becomes corrupt or fills as completely blank. Here is an example of a corrupted image (colored right-side of main map)

I'm using the OpenLayers Statmen plugin for the basemap. I've tried a few others and they either fail as well or do not allow the map tile to be printed at the scale. If it clear cache and reboot it will work for about 7-10 images and then fail again.


texts: QGIS generates multiple WMS requests in single image WMS mode

texts:
I am testing my WMS service with QGIS and have noticed that even when the layer is configured in "non tile" mode (tile size is not specified) QGIS will generate multiple requests to my WMS service. This will only happen when the map area is large (QGIS running on large resolution monitor).

Is there a way for QGIS to only generate one request WMS service (with image pixel size request equal to exactly currently viewed map size)?


texts:
I recently tried to update my QGIS version with the installation of Osgeo4w libraries. I was told that the installation was supposed to take something like 20 gb space on my computer.

While downloading the libraries and not even crossing the half way line, a 20 gb library space has already been downloaded, and it didn't even showed sign of finishing.

I tried to find the location of the libraries on my computer, but I couldn't find it. The path was supposed to be a temporary folder called "C:\Users\CARMEL~1\AppData\Local\Temp", but it just didn't exist.

Any body know where I can find those libraries? Is there a new path where it goes as default after finishing its download?


texts: How to perform a cascaded combine from python script in QGIS?

texts:
I'm working with a large dataset, and I need to be able to combine up to 400,000 polygons.  I have been using the QgsGeometry::combine(QgsGeometry* geometry) function and simply looping over all geometries in an interative fashion like this:

union_geoms = geometries[0]
for geometry in geometries:
  union_geoms = union_geoms.combine(geometry)

This works fine for most of my data, (the largest list of polygons I have gotten it to complete is ~30,000.)  But I have two sets of ~200,000 and one of ~400,000 and I can't get this to complete in a timely fashion (I have let it run for days to no avail).  So I started looking for a better method to combine these geometries, and I found geos::operation::geounion::CascadedPolygonUnion  This algorithm, described here, is well suited to my application as the polygons of interest contain a lot of adjacency.

I found a couple references to implementations of this in Qgis, here in 2010 and here in 2013, but it doesn't appear to have made it into the API (unless I'm looking in the wrong place!)

So my question is how would I go about implementing this type of algorithm, or more simply how can I use the existing implementation in GEOS from my own (python) scripts?  From the 2013 link,

  you can use QgsGeometry::asGeos() in the plugin instead of wkt-export / wkt-import. For cascaded union, there is the method GEOSUnaryUnion(const GEOSGeometry* g1) in the c-api. The argument seems to be a geometry collection, which can be created with GEOSGeom_createCollection.

But there doesn't seem to by python bindings for the asGeos() method.

Anyone have suggestions as to how to proceed?


texts: QGIS Conditional statements in field calculator

texts:
I have attempted to write an expression for hours now by looking through these forums.

I want to write if this column equals Blank, than this column equals blank.
      Ex) If defect_typ = A_FC THEN Defect_Nam = Fatigue Cracking

I can get a return from multiple codes, but here is an example of one I did:

CASE WHEN defect_typ LIKE '%A_FC%'
       THEN 'Fatigue Cracking' ELSE 'other'

No matter what code I use it fills my entire defect_nam column with Fatigue Cracking or other, and if I try to do multiple, (which I need to create multiple conditional rules) it will fill the column with a mash up of the returns (fatigueCrackingBlockcrackingPatching) but I need to figure out how to do one before I can move on.

I  am an ArcMap Vet and to get the result I want in Esri I would use the following script

def Defect(defect_types_code):
    if (defect_types_code == "A_LC"):
      return "Linear Cracking"
    elif (defect_types_code == "M_PS"):
        return "Parking Bumper Damage"
    elif (defect_types_code == "M_RB"):
        return "Rebuild/Replace Catch Basin"
    elif (defect_types_code == "M_RD"):
        return "Ramp Damage
    else:
      return "N/A"

How do I write this in QGIS?

I was linked to another question with answers to help, "Elseif Conditional Statement in QGIS Field Calculator,"but the outcome of this code filled my defect_nam with all Nulls.


texts: Errors opening Geospatial PDF in QGIS 3.8.3

texts:
I'm a complete QGIS newbie so bear with me if I'm missing something obvious...

I'm unable to open GeoPDF's that I've downloaded from USGS's The National Map (TNM) site. This is the error that I get in the QGIS Log:

2019-09-24T12:44:03     CRITICAL    Invalid Layer : GDAL provider Cannot open GDAL dataset /Users/jeff/Downloads/IL_Englewood_20120827_TM_geo.pdf:
             `/Users/jeff/Downloads/IL_Englewood_20120827_TM_geo.pdf' not recognized as a supported file format.
             Raster layer Provider is not valid (provider: gdal, URI: /Users/jeff/Downloads/IL_Englewood_20120827_TM_geo.pdf

When I run gdalinfo on the file I get this:

gdalinfo ~/Downloads/IL_Englewood_20120827_TM_geo.pdf
ERROR 4: `/Users/jeff/Downloads/IL_Englewood_20120827_TM_geo.pdf' not recognized as a supported file format.
gdalinfo failed - unable to open '/Users/jeff/Downloads/IL_Englewood_20120827_TM_geo.pdf'.

From everything I've seen online GeoPDFs should be fully supported in QGIS 3.8.3 and gdal 2.4.1. I can also open the file in the standard Apple Preview app without trouble.

Here's my QGIS "About" output if that's helpful:


texts: Creating an image collection from latest computed median images in Google Earth Engine

texts:
I want to create an image collection where all the median images can get stored so that later on it can be exported to drive.

Link to my code: https://code.earthengine.google.com/4dc5685ed79044086587d811c4ee584f

//cloud mask along with function to add image as layer in the collection
function addImage(img)
{
var id = img.id;
var image = ee.Image(img.id);
Map.addLayer(image);

var cloudBitMask = ee.Number(2).pow(10).int();
var cirrusBitMask = ee.Number(2).pow(11).int();
var qa = image.select('QA60');
Map.addLayer(qa.bitwiseAnd(cloudBitMask).neq(0), {}, 'clouds');
Map.addLayer(qa.bitwiseAnd(cirrusBitMask).neq(0), {}, 'cirrus');

function maskS2clouds(image)
{
var qa = image.select('QA60');
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
         qa.bitwiseAnd(cirrusBitMask).eq(0));
return image.updateMask(mask);
}

var cloudMasked = sentinel.filterBounds(agra).map(maskS2clouds);
var median = cloudMasked.median();
Map.addLayer(median
, {bands: ['B4', 'B3', 'B2'], max: 2000}, 'median');
print(median);
}


texts: Converting Raster Orthography from MrSID to GeoTIFF Creates Colourful Artifacts

texts:
I'm using geoserver 2.9.0 to display our orthography. I've managed to get the MrSID file displaying properly however it's very slow. From what i've read to increase the performance I should convert to a GeoTIFF that is tiled and has overviews. My first step is to decode the MrSID file and save it as a geoTIFF.
I do that with the following command:

mrsidgeodecode -i Mosaic_10cm_3TM.sid -o Mosaic_10cm_3TM.tif -of tifg -progress timer

*mrsidgeodecode is from lizardtech: https://www.lizardtech.com/gis-tools/tools-and-utilities#hs_cos_wrapper_widget_1460663711843

Here are some samples from the imagery served by geoserver over WMS:

MrSID Building (Original):

GeoTIFF Building (Converted):

MrSID RV Dealer (Original):
http://i.stack.imgur.com/JIbuu.jpg

GeoTIFF RV Dealer (Converted):
http://i.stack.imgur.com/C54es.jpg

Following the decode process I use the following GDAL commands to build the raster into a file that is much more performant:

REM -- Reproject GeoTif from EPSG:3780 to EPSG:3776 (http://www.gdal.org/gdalwarp.html)
gdalwarp -s_srs EPSG:3780 -t_srs EPSG:3776 -r bilinear -multi -wm 4096 -of GTiff Mosaic_10cm_3TM.tif Mosaic_10cm_3TM_EPSG3776.tif

REM -- Tile GeoTif (http://www.gdal.org/gdal_translate.html)
gdal_translate -of GTiff -co TILED=YES -co BIGTIFF=YES Mosaic_10cm_3TM_EPSG3776.tif Mosaic_10cm_3TM_EPSG3776_tiled.tif

REM -- Add Geotif Overviews (http://www.gdal.org/gdaladdo.html)
gdaladdo --config BIGTIFF_OVERVIEW YES --config GDAL_TIFF_OVR_BLOCKSIZE 256 -r gauss Mosaic_10cm_3TM_EPSG3776_tiled.tif 2 4 8 16 32 64

All of these GDAL commands work and create an image that is much more performant but because the decode process creates the colourful artifacts they remain present throughout.

Oddly enough however when I create an ImagePyramid with the following command (using the source that contains the colourful artifacts *created by mrsiddecode):

gdal_retile -v -r bilinear -levels 6 -ps 2048 2048 -targetDir pyramid Mosaic_10cm_3TM.tif

I get an image that is nearly identical to the original MrSID file:
http://i.stack.imgur.com/4aB3T.jpg

Any help removing the colourful artifacts or converting the data in a different way as to not produce the artifacts would be great!


texts:
I am currently editing the Configurable Map Viewer so that it can display my map service from local ArcGIS Server.

I can't find which HTML element is the container for the map service. I can not find "DIV id = map".

Does a map to be displayed always contained in DIV, or there is another method to display the map on the page?


texts: Replacing raster values in one raster with another raster in QGIS ... how to?

texts:
I have a raster (DEM) and want to replace a few cell's values with the values from another raster in QGIS. I would have done this with the "over" tool in the raster calculator in ArcGIS.

Specifics: To potentially restore a mire, I want to model how the water flows if dams were to be put into place across some ditches. So, I have the original DEM and would like to add the dams to it (to then do an flow accumulation analysis). I cannot just add a few meters to the cell's values, because the dams are supposed to go across the ditches (so 90 degrees to the ditches) ... so, a constant would not result in an even dam, but a v-shaped one.

Is somebody able to help?


texts:
I found there are a number of RDMS have the spatial capability e.g. mysql, postgreSQL

I considering which I should go forward!

I would like to store the lat lon coord from the tracker to the db and make a query to find the nearest place or compare with the defined location.

Could anyone give me an idea on this topic?


texts: How to Create dojox.form.CheckedMultiSelect dropdown?

texts:
I'm trying to create a drop down multi select checkbox using dojo JS and html but just given the drop down only.
Code:
<body>
<div>
 <select data-dojo-type="dojox/form/CheckedMultiSelect" data-dojo-props="dropDown:true,labelText:'Countries',multiple:true, required:false" name="country">
            <option value="THA">Thailand</option>
            <option value="VEN">Venezuela</option>
            <option value="KWT">Kuwait</option>
            <option value="FRA">France</option>
            <option value="CHL">Chile</option>
        </select>
</div>
</body>
<script>
define([
  'dojo/_base/declare',
  'jimu/BaseWidgetSetting',
  "dojo/dom",
  "dijit/form/TextBox",
  "dojo/dom-construct",
 ,"dojo/_base/lang","dojo/on","dojo/ready","dojo/dom-style",'jimu/dijit/Message',"esri/request","dojo/promise/all","dojox/form/CheckedMultiSelect","dojo/parser",
  "dojo/domReady!"
],
function(declare, BaseWidgetSetting,dom,TextBox,domConstruct,FeatureLayer,lang,on,ready,domstyle,Message,esriRequest,all,CheckedMultiSelect,parser) {}
</script>
//@css....................................................
@import "//js.arcgis.com/3.22/dojox/form/resources/CheckedMultiSelect.css";
.dojoxCheckedMultiSelectCheckBoxInput {
    opacity:1
}
.dojoxCheckedMultiSelectSelect{
  width: 240px;
}
.dojoxCheckedMultiSelectWrapper {
    width: 200px;
}


texts: Geopandas .difference between two different shapes produces empty difference

texts:
I have been puzzled by this, and I feel like it's a geopandas bug.
base = gpd.GeoSeries(dist_to_remove).plot(edgecolor='black', color='yellow')
gpd.GeoSeries(prec_shape).plot(ax = base, color='red', edgecolor='black')

So this should be an easy difference to compute. However, upon doing
dist_to_remove.difference(prec_shape)

I get the above (there is a tiny tiny square in the lower left). The area has shrunk by an enormous magnitude. I am really confused as to what is happening. The odd part is, if I change the buffer on prec_shape, the difference occurs as expected. I'm really confused as to what's going on right now, is this a shapely/geopandas bug?


texts:
I used the instructions at https://pythongisandstuff.wordpress.com/2011/07/07/installing-gdal-and-ogr-for-python-on-windows/

I downloaded the following files:

GDAL-1.11.1.win-amd64-py2.7.msi
gdal-111-1800-x64-core.msi

Then I followed the above instructions.  When I type

from osgeo import gdal

at the Python IDLE command line, everything appears to be OK. Then I type in

ds = gdal.Open( 'C:/testfile.DT1')

and again, it appears to be OK. However, when I type in

gt = ds.GetGeoTransform()

I get the following error message:

Traceback (most recent call last):
  File "<pyshell#3>", line 1, in <module>
    gt = ds.GetGeoTransform()
AttributeError: 'NoneType' object has no attribute 'GetGeoTransform'

I have checked the DTED file and there is nothing wrong with it.  I also tried the above with a tif and got the same result. I can get the GDAL executables to run, but not the GDAL Python scripts. Is there something that I should check regarding my install of GDAL?


texts: Adding slope values to attribute table (QGIS 3)

texts:
I have a DEM raster and some points (one .shp for all points) plotted on it.

I would like to know if there is a way to add the slope value of the pixel that each point is located in the attribute table of the point's shp.

Solutions using Python are also welcome.


texts: Extracting NDVI by regions defined by certain column using Google Earth Engine

texts:
I am trying to extract NDVI values by regions of a certain polygon shapefile. I need, as output a graph (or a table) where the regions are defined by a certain column, called 'FID_incend'. I need exactly this column because those are unique values that I will use afterwards to process other maps.

My code is:

var incendis1987 = ee.FeatureCollection("users/oriolbaena/incendis1987_identityLitology_87929702_s_0p9HA_simplified");

Map.centerObject(incendis1987, 12);
Map.addLayer (incendis1987,{color: 'FF0000'})

var addNDVI = function (image){
  var ndvi= image.normalizedDifference(['B4', 'B3']).rename('NDVI');
  return image.addBands(ndvi)
};

var dataset = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR")
    .filterBounds(incendis1987)
    .filterDate('1984-01-01', '1995-12-31')
    .map(addNDVI).select('NDVI');

var graph = ui.Chart.image.seriesByRegion({
imageCollection: dataset.select('NDVI'),
regions: incendis1987,
reducer: ee.Reducer.mean(),
scale: 30,
seriesProperty: FID_incend,
})
  .setChartType('ScatterChart')
  .setOptions({title: 'NDVI Incendis1987 dates1984-1995', vAxis: {title:'NDVI'},
             lineWidth: 1,pointSize: 2,});
print(graph)
print ('I love JavaScript (L)')

As you can see, attibutes of my table include the column 'FID_incend':

The shp can be found here: https://code.earthengine.google.com/?asset=users/oriolbaena/incendis1987_identityLitology_87929702_s_0p9HA_simplified


texts: GetFeatureInfo of a specific GeoServer layer by selecting it from a set of GeoServer layergroup

texts:
I have multiple Geoserver layers, suppose 2 here (Rainfall and Maximum Temperature)
and I want to extract the feature of that layer which I check/switch in the sidebar,
But for now I just get the feature of rainfall,
as the URL source (var url) that I am passing in .JS is of rainfall only,
var url = india_dist_rainfall_layer_source.getFeatureInfoUrl

how to make this URL source dynamic in GetFeatureInfo function of Openlayers 6, to get the info of that specific layer which we select from sidebar?

Weather Parameters

                <input type="radio" name="WeatherParameterLayerRadioButton" value="Rainfall" checked> Rainfall

                <input type="radio" name="WeatherParameterLayerRadioButton" value="Maximum_Temperature"> Maximum Temperature

window.onload = init;

function init(){

    // 1.1 Rainfall
    var india_dist_rainfall_layer_source = new ol.source.TileWMS({
        url: "http://localhost:8080/geoserver/agrodss/wms",
        params: {"LAYERS":" agrodss:Rainfall (mm)", "tiled": true},
        serverType: "geoserver",
    })
    var india_dist_rainfall_layer = new ol.layer.Tile({
        source: india_dist_rainfall_layer_source,
        opacity: 0.0,
        visible:true,
        title: "Rainfall"
    })

    // 1.2 Maximum Temperature
    var india_dist_maximum_temperature_layer_source = new ol.source.TileWMS({
        url: "http://localhost:8080/geoserver/agrodss/wms",
        params: {"LAYERS":"agrodss:Maximum Temperature (°C)", "tiled": true},
        serverType: "geoserver",
    })
    var india_dist_maximum_temperature_layer = new ol.layer.Tile({
        source: india_dist_maximum_temperature_layer_source,
        opacity: 0.0,
        visible:true,
        title: "Maximum_Temperature"
    })

    var myview = new ol.View({
        center: ol.proj.fromLonLat([80, 22]),
        zoom: 3,
        maxZoom: 9,
        minZoom: 2,
      });

    var map = new ol.Map({
        target: 'js-map',
        view: myview,
    });

    var weatherParameterLayerGroup = new ol.layer.Group({
      layers: [
        india_dist_rainfall_layer,
        india_dist_maximum_temperature_layer
      ]
    })

    map.addLayer(weatherParameterLayerGroup);
    var weatherParameterLayerElements = document.querySelectorAll('.sidebar > input[type=radio]');
      for (let weatherParameterLayerElement of weatherParameterLayerElements){
        weatherParameterLayerElement.addEventListener('change', function(){
          let weatherParameterLayerElementValue = this.value;
          weatherParameterLayerGroup.getLayers().forEach(function(element, index, array){
            let weatherParameterLayerTitle = element.get('title');
            element.setVisible(weatherParameterLayerTitle === weatherParameterLayerElementValue);
          })
        })
      }

    map.on('singleclick', function (evt) {
        document.getElementById('rainfall_info').innerHTML = '';
        var myviewResolution = myview.getResolution();
        var url = india_dist_rainfall_layer_source.getFeatureInfoUrl(
        evt.coordinate,
        myviewResolution,
        {'INFO_FORMAT': 'text/html',
        },
        );
        if (url) {
        fetch(url)
        .then(function (response) { return response.text(); })
        .then(function (html) {
        document.getElementById('rainfall_info').innerHTML = html;
          });
        }
    });
}


texts:
I have an issue for the result of QGIS Zonal statistics  mean values.
I expect the value for the polygon marked red (poquoson) is about 26.6 but it seems that it is averaged value over one right shifted cells.
Values showing in small black font are mean values for the each polygon
according to zonal statistics.
Values showing in blue are raster values in us.tmax_nohads_ll_20140525_float.tiff

Here is "identify results"

I made zonal statistics using files with
tl_2013_us_county.shp   &    us.tmax_nohads_ll_20140525_float.tiff
on http://www.qgistutorials.com/en/docs/sampling_raster_data.html

1st thing I did was to change CRS for tl_2013_us_county.shp  to EPSG4326.
I believe that CRS of us.tmax_nohads_ll_20140525_float.tiff  and
tl_2013_us_county.shp became the same in such a way.

2nd I executed zonal statistics

Platform is QGIS2.18.13 32bit
On Windows 10 Pro 64bit Intel Core i5-6500 CPU

Is there any way to get values what I expect?

For above case, expected value is over 26.06 instead of 25.08.


texts: R raster extract() and keep the original shapefile attributes

texts:
I use the R function extract() in order to extract the values from a raster dataset. The shapefile that I use with the extract() function has some columns/attributes -for example the unique ID and a characteristic of each polygon- that I need to remain within the extracted dataframe.
The readOGR function returns a SpatialPolygonsDataframe that has a "data" dataframe and a "polygons" list.
I would need the extracted dataframe with extract(raster, shapefile, df=TRUE) to associate the extracted values with the shapefile "data" attributes, also to be sure of what is really what (maybe someone can explain -if it is relevant- in order to avoid possible confusion, about the attribute from readOGR named "plotOrder").

A sample code:

require(sp)
require(raster)

setwd('FILEPATH')
shapefilepath <- 'SHAPEFILEPATH'
rasterpath <-'RASTERFILEPATH'
shp <- readOGR(dsn = ".", layer = "SHAPEFILENAME_withoutDOTshp")
pts <- slot(shp, "data")
img <- brick(rasterpath)
rastervalues <- extract(img, shp, cellnumbers=TRUE, df=TRUE)

So, this returns a dataframe that does not include the associated "data" from the SpatialPolygonsDataframe from readOGR().
Any way to overcome this?


texts: I've installed QGIS 2.16 and have noticed that the Merge (Vector Menu>Data Management) doesn't appear to have a Folder option anymore. This was the option that allowed you to process/merge a complete folder without having to load layers into your project/workspace first. I've had a look at the changelog for QGIS 2.16 - although fTools have been replaced by Processing tools I can't see why the extra functionality has been removed.


texts: How to symbolize features with NULL values in graduated symbology?

texts:
I have a polygon feature dataset, and an attribute from a different table that I join (one-to-one) within QGIS in order to symbolise the attributes as a choropleth map. However, not all polygon fields have a matching field in the table of numerical attributes, so there are some null values when producing the graduated symbology.

From my research, the most common piece of advice in this situation is to include a copy of the polygon dataset (or some other background) that defines a default symbology. This works because the null values are not classified in the graduated symbology, so one can see "underneath" them. I have attached an image of exactly this. The dark grey features do not exist in the join table (mb_percentile_isochrones_all), but do exist in the boundary polygons table (mb2013_wgtn). So I need two instances of the mb2013_wgtn table in order to show the "no data" features.

However, this does not seem elegant to me. Much more intuitive would be to define a null value symbol. Perhaps this could be perfectly transparent to be consistent with what currently exists, or perhaps it would be some kind of muted grey—whatever the user wants. The point is, at present one needs two different layers in the contents in order to handle the symbology of null values. This means that to change the symbology of all your features at the same time (say, if you want to increase the width of all borders), this has to be handled twice: once in the properties for the feature with the graduated symbology, and once for the "background" layer that handles the null symbol.

Is it possible to define a "null" value symbol without using a "copy" of the same layer used for the graduated symbology, in QGIS (2.6.1)?


texts: Global name 'dataobjects' is not defined in QGIS?

texts:
I have this small problem with QGIS 2.18.4 Las Palmas. I just downloaded it and tried to merge vector layers, but i've got error: "global name 'dataobjects' is not defined".
I didn't change anything that in my opinion could cause this problem. Although I wanted to use lastools via QGIS, so I connected it in previous version (as far as I remember it was QGIS 2.14) and due to some errors I had to remove some folder paths to Grass, but everything worked just fine, so it's probably not the cause.
I tried to merge shapes by the plugin, but it didn't work either.
Any ideas how i can fix this?


texts: How to use rio command line from rasterio?

texts:
I want to use rasterio and rio-color plugin.

I used pip for installing rasterio and rio-color and afterwards when I want to run "rio", I see this:

 Traceback (most recent call last):
  File "/usr/local/bin/rio", line 7, in <module>
    from rio.manage import main
  File "/usr/local/lib/python2.7/dist-packages/rio/__init__.py", line 5, in <module>
    from .app import create_app
  File "/usr/local/lib/python2.7/dist-packages/rio/app.py", line 11, in <module>
    from rio.setup import configure_app
  File "/usr/local/lib/python2.7/dist-packages/rio/setup.py", line 10, in <module>
    from .core import db
  File "/usr/local/lib/python2.7/dist-packages/rio/core.py", line 14, in <module>
    from flask_user import UserManager
  File "/usr/local/lib/python2.7/dist-packages/flask_user/__init__.py", line 9, in <module>
    from flask_login import LoginManager, UserMixin as LoginUserMixin, make_secure_token
ImportError: cannot import name make_secure_token

What should I do?


texts: Raster conditional tool unexpected results. What's wrong?

texts:
I have a DEM raster (High: 53441, Low: -3444, unit: cms) with negative values, and I want to make all values 100 or less to be 100.

I've tried Raster Calculator, and also python. I've got the same unexpected result: (High: 22881, Low: 8842)

This is my code in python:

inRaster  = Raster('ns_dem_int')
outRaster = Con(inRaster, 100, inRaster, "Value <= 100")
outRaster.save("DEM14")

I've run the code -and the Raster Calculator- several time (with some variations) but I've always end up with those numbers as results
Any ideas on what could be wrong here?


texts: Creating a sidewalk map for a walkability survey

texts:
If I wanted to do my own walkability survey of my city, what tools could I use? I have map data from OpenStreetMap, but I have to consider that streets have two sides (which aren't always accessible from the other side of the street). I guess I need to make a sidewalk map.

I'd like to use path-finding algorithms against the data to provide walking directions that avoid sections of streets with poor walkability. Some blocks have sidewalks, some don't, and some streets have a sidewalk on one side but are too busy to cross. The directions would cross the highway at places where it's easier--very easy at the crosswalk, hard at several other streets, impossible at a few streets.

I know that weighted path-finding algorithms exist. I'm just wondering how to create this data. I've toyed with taking a map image, putting it in Inkscape, turning on grid lines, and overlaying paths along the sidewalks. With the grid lines, the points of paths connect for path-finding. Path color indicates difficulty of walking. Aside from the path-finding, it'd be a good map to show walkability in town.

Anyone have any ideas?


texts: Publishing GeoTiffs with Geoserver - ColorPalette with Alpha

texts:
I want to publish rasterdata (Geotiffs) with GeoServer. The resulting image from WMS is shown in Grayscale, but actually it is an RGB-Image.
gdalinfo:

P:\OSGeo4W\bin>P:\OSGeo4W\bin\gdalinfo  Q:\data\tif\pic.tif
Driver: GTiff/GeoTIFF
Files: Q:\data\tif\pic.tif
Size is 4523, 4522
Coordinate System is:
LOCAL_CS["DHDN / 3-degree Gauss-Kruger zone 3",
    GEOGCS["DHDN",
        DATUM["unknown",
            SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    AUTHORITY["EPSG","31467"],
    UNIT["metre",1]]
Origin = (3532474.929505631800000,5367000.280343170300000)
Pixel Size = (0.333485342299112,-0.333485342299112)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2012/05/08 11:38:04
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=WGEO 5.0 WGEO BASIS
  TIFFTAG_XRESOLUTION=200
  TIFFTAG_YRESOLUTION=200
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 3532474.930, 5367000.280)
Lower Left  ( 3532474.930, 5365492.260)
Upper Right ( 3533983.284, 5367000.280)
Lower Right ( 3533983.284, 5365492.260)
Center      ( 3533229.107, 5366246.270)
Band 1 Block=256x256 Type=Byte, ColorInterp=Palette
  Overviews: 2262x2261, 1131x1131, 566x566
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 2262x2261, 1131x1131, 566x566
  Color Table (RGB with 256 entries)
    0: 0,0,0,255
    1: 128,0,0,255
    2: 0,128,0,255
    3: 128,128,0,255
    4: 0,0,128,255
   ...
  253: 255,0,255,255
  254: 0,255,255,255
  255: 255,255,255,255
Band 2 Block=256x256 Type=Byte, ColorInterp=Alpha
  Overviews: 2262x2261, 1131x1131, 566x566

I think the problem is the Alpha-Channel. Without this Channel, the resulting image/colors are correct.
I'm using the default style "raster" from Geoserver.


texts: Determine if a polygon intersects itself in OpenLayers

texts:
I let the users draw polygons and then I try to save them out to the SQL Server 2008 SQLGeography field, but it is giving me an error:

  Error 24200: The specified input does not represent a valid geography instance.

I think the issue is because the polygon intersects itself and SQL doesn't accept this behavior. Is it possible to do a check in OpenLayers to see if the WKT from the polygon would be accepted by the SQLGeography.STGeomFromText function?


texts:
I have some long linear polygons that I'd rather store as lines but I'm not sure how best to do this. Ideally I'd have a line going through the centre of the polygon lengthways. I have two ideas of how to do this but I'm not sure how to implement them or to what extent they're possible. Perhaps there's already a tool out there that does this?

option one:
1. convert polygon to points at vertices
2. identify the two points furthest apart in each feature. I can do this in point distance but I'm not sure how to automate it
3. mask the area outside the polygon
4. run a cost path between the two points

option two:
1. convert polygon to points at vertices
2. Create theissen polygons around the points
3.(somehow) save the line along the middle between the theissen polygons?

Does anyone know how I can make either of these methods work or if there is an easier method. My Python is weak so I was going to try it in model builder but if I needs Python I can give it a go


texts: FME - How to extract the rightmost folder in a path

texts:
I have an FME workbench which reads all of the TAB files in a folder and all of the subfolders and writes the paths of those TAB files as well as their type (NATIVE, RASTER, etc) to an Excel file.

The folder is a published parameter.  What I would like to do is name the Excel file the same as the rightmost folder (I'm not sure what you'd call this) in the parameter.

For example, when the workbench is run, the user selects this folder:

C:\Temp\ABC\DEF\GHI

The excel file would be GHI.xlsx

or if the user selects this folder:

C:\Folder\ABC\DEF123

The excel file would be DEF123.xlsx

I have tried many different ways to get this folder name.  I tried creating a new parameter that builds the name but no luck.  I've also tried setting the parameter as an attribute using the parameterFetcher then using the SubstringExtractor.  However, I can't seem to get the start and end index set properly.  This is what I've got in my substringExtractor:


texts: Clipping points using buffer in Arcmap takes FOREVER

texts: I'm clipping points (crimes) by using a buffer around a shapefile of city parks. There are quite a few records involved (2 million+), but it is still taking forever - as of now, the clip has been running for an hour and ten minutes. I'm also using a subset of the parks data, so I'm not clipping crimes with every polygon in the shapefile.  My questions: How does Arcmap go about it's clipping? Does it have to check every point, even those that are nowhere near the polygons? Is there a better method? I have QGIS at home, but I'm not sure if it's as fast as Arcmap. Thanks, Ward


texts: Copying similar files from geodatabases into single geodatabase using ArcPy?

texts:
I'm trying to copy files form multiple geodatabses into a single database. I'm looking for two subclasses. Here's the code so far, it isolates the feature classes but when I try to run the search cursor it only runs it on the last feature class that the arcpy.da.Walk throws out. Any ideas on how I can fix this?

import arcpy, os, string, sys
workspace = r"C:\Users\A\Desktop\A"
outWorkspace = "C:\Users\A\Desktop\A"
arcpy.env.workspace= workspace

for dirpath, dirnames, filenames in arcpy.da.Walk(workspace, datatype="FeatureClass"):
    for filename in filenames:
       arcpy.env.workspace=filename
       desc = arcpy.Describe(dirpath)
       featureTest=(os.path.join(dirpath,filename))
       print filename

fc= filename
fields = ["F_CODE"]

with arcpy.da.SearchCursor(fc, fields) as cursor:
    for row in cursor:
        print (row)


texts: Should I query the tab/shp file or SQL (mapserver or geoserver)?

texts:
I have tab files (or shp) with information like residents , buildings etc...
I need to be able to click a point in web map and query all data that fall in certain radius.
I didn't decide yet which WMS I'll use because my web apps are .net and although geoserver seems more popular , I have zero experience with java.
Let's say I chose a WMS. Should I use it's sdk (if it has any , currently using mapxtreme to query map data on server side) to query files directly or should I import all to SQL 2008 and work with it (not sure how to search within polygon or in radius yet)
Tips , opinions will be welcome


texts: ArcGIS Server - Can python GP services returns json data to client through REST?

texts:
Can python GP services returns json data to client through REST?

Imagine I have a python GP service that creates buffers around points and then merge the polygons together. I would like to send the output in esri json format as the response of the REST request made by the javascript client.

The client could use ajax to make the REST request. How would I make the GP to return the json to the client when job is finished?


texts: Creating polygons from centroids in Python for an interactive map?

texts:
I have a geopandas dataframe where each record is a point. I am making an application using Flask and Python and would like to display an interactive map in HTML where the points would be transformed into a raster or into polygons from centroids. Each pixel would have an onclick event. I used geopandas simplify and dissolve and buffer (separately) to create polygons from the points with ok-looking results but then I can only display the map as a static image in HTML. I tried Bokeh but the best I can get is square points on a map which do not look good because I want them to be touching as a continuous structure but do not know the scale in advance so cannot set the right size of the point.
Ultimately, I would like it to be like a grid in this tutorial (but they use polygons): https://automating-gis-processes.github.io/2016/Lesson5-interactive-map-bokeh.html

Would anyone have an idea how to go about this?


texts: C# Library for consuming Geoserver REST?

texts:
I am currently experimenting with Geoserver and it's capabilities within C#. One thing I would like to do is query the server for basic information (existing workspaces, layers, layer attributes, etc...)

I noticed that Geoserver does have a REST implementation, but is there a C# implementation for consuming this service? I am sure a library could easily be written, but I would hate to re-invent the wheel if one already exists. Thank you!

Edit

This can be done using C# WebRequests. However, what I would have created is a C# library for consuming the Geoserver REST implementation. It seems like something many people would use. Another way to ask the question would be, "Does a library full of C# WebRequests exist to consume Geoserver's REST implementation?"

It is much easier if a library already exists with classes such as these to consume the webservice, than it is to create an entire library. For Example:

var geoserver = new Geoserver("http://localhost:8080/geoserver/web");
var workSpaces = geoserver.Workspaces;
var layer = geoserver.Workspaces[0].Layers["Foo"];

The above code is the type of implementation I would really enjoy. If I was to create a library that consumed the Geoserver REST service. The above is what I would create for myself to use. If it wasn't created already, I would likely make an opensource project that others could use as well. So that in the future, if someone is wanting a library to easily consume the Geoserver service, they can contribute and use the one that already exists! :)


texts: arcpy path of ArcGIS mapservice layer throws error

texts:
I have a few basemaps published to ArcGIS Server. I have a script that I am running which removes these basemap map services. What I need to do is find the source REST endpoints of these mapservices to further narrow down my script. I have just added as per @Vince suggestion the isServiceLayer as an if statement and tried printing out the serviceProperties but the script throws an error on certain layers

sample UPDATED code...

f = path
mxd = arcpy.mapping.MapDocument(r'{}'.format(f))
for df in arcpy.mapping.ListDataFrames(mxd):
    for x in arcpy.mapping.ListLayers(mxd,data_frame=df):
        if x.isGroupLayer:
            name = str(x)
            if name.startswith("Basemaps"):
                if x.isServiceLayer:
                    print x.serviceProperties
                    arcpy.mapping.RemoveLayer(df,x)

ERROR

_Basemaps_hcgis
_Basemaps_hcgis\ Region
{u'UserName': u'arcgis', u'Supportsrsa': u'-1', u'ServiceType': u'MapServer', u'Connectionfile': u'C:\\Users\\\\AppData\\Roaming\\ESRI\\Desktop10.3\\ArcCatalog\\arcgis on SERVER (admin).ags', u'Connectionmode': u'1', u'URL': u'http://SERVER', u'Servertype': u'1', u'Restversion': u'10.3', u'Rsaurl': u'http://SERVER', u'Server': u'', u'Adminurl': u'http://SERVER', u'Resturl': u'http://SERVER', u'Soapurl': u'http://SERVER', u'Handshakefinished': u'-1', u'Anonymous': u'0', u'Usedefaultstagingfolder': u'-1', u'Admintokenurl': u'http://SERVER', u'Stagingfolder': u'C:\\Users\\\\AppData\\Local\\Temp\\arc434C\\Staging'}
_Basemaps_hcgis\ Region\>570000
Traceback (most recent call last):
  File "C:\Users\Desktop\blah.py", line 11, in <module>
    print x.serviceProperties
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\arcobjects\_base.py", line 78, in _get
    (attr_name, self.__class__.__name__))
NameError: The attribute 'serviceProperties' is not supported on this instance of Layer.

UPDATE from @Alex suggestions:

my mapserver basemaps are in a group layer named Basemaps and I only want to remove those basemaps... the problem is I am running this script through network drives so some of the old MXDS(before my time) have two group layers called Basemaps group layer, 1 from the old server and 1 from the new one. I DO NOT want to remove the new **Basemaps ** group layer and the only way to accomplish this is to open up each mapservice and check which server it is from... here is a picture of how it is set up


texts: How to create line from parallel line objects

texts:
I have a two tables of lines in SQL Server 2016 database. Each one has a Spatial field with type geometry.
My job is to find lines from table B which are parallel line objects to lines in table A and make line object from them.
The thing is I'm new at spatial SQL and I really can't figure out how to do that.
My tables look's like this:

My main line look's like this, table A:

And example of parallel line objects, table B:

Here is desired result which I'm looking for:

I've tried something with ShortestLineTo and delimited lines from table B to start and end points but without any success. Literally any help is welcomed.


texts: Which horizontal unit is used on the table of the Profile Tool?

texts:
Does anyone knows what horizontal unit is used on the table of the "Profile Tool".
Vertical is in meter but I don't know about the horizontal axis.

Knowing the length of my profile, I tried to divide that length by the number of extracted elevation but I noticed that the gap between each elevation data seems to change.

At first I have one elevation data each 0.00453 for a hundreds of data. Then its one data per 0.0060 then 0.0105 etc etc.

Would that mean that elevation data are not extracted at equally spaced points ?

Thank you very much in advance !

Note: I am using as raster a .tif DEM and a line built in shapefile.


texts: Finding Land Surface Temperature in across 4 continuous images in ArcGIS 10

texts:
I'm trying to find Land Surface Temperature of 4 continuous Landsat8 images in ArcGIS. I was able to calculate them as individual panels (using these two videos as guidelines: https://www.youtube.com/watch?v=uDQo2a5e7dM ;
https://www.youtube.com/watch?v=GwwXhMKo-yM), but that created nodata values that I couldn't eliminate.

My alternate approach was to create a mosaic of the images using composites of each image, but this leaves out bands 10 and 11, which I need to do the LSE calculations.  The concern I have with this is that all four images have different thermal band values/properties; if these bands don't make it into the mosaic in the way the images will, how can the calculation be done?

How can I do this in ArcGIS Desktop?


texts: Nodata value making holes in merged GeoTIFFs in QGIS?

texts:
I am trying to remove black borders around the edge of a GeoTIFF. The Black border around the GeoTIFF has an RGB value of (0,0,0). I have tried assigning 0 as a "No Data" value to the tiff, but QGIS then makes all pixels which contain a 0 such as (0,54,76) into "No Data" pixels. This leaves holes in the map.

Is there a way of only classifying a "No Data" pixel if 0 is found across all 3 RGB bands on a pixel(pure black pixel)?

I am using the latest version of QGIS 3.8 Zanzibar. Here are 2 screenshots to clarify the situation.

Image of file for merging