texts: Selecting ArcSDE polygon by point in ArcGIS Desktop using ArcPy?

texts:
I keep thinking that I must be missing something, but there does not seem to be a tool in ArcGIS 10 to select features (in particular polygons) from a layer at a point (X,Y) location via ArcPy.  The parameters for such a tool would just be a layer name and an XY location.

At the moment I workaround this by creating a point featureclass containing the point and performing a SelectLayerByLocation on it.  However, when the polygon feature class is in Oracle (accessed via ArcSDE 9.x) and contains 3.5 million polygons the time taken to make the selection can be more than 5 mins when I think a second or two (with less code) would be more appropriate.  The feature class has a spatial index and I've tried using arcpy.env.extent (which SelectLayerByLocation appears to ignore) to restrict the geographic area accessed but the performance remains very poor.

Is there a quicker way to do this using ArcGIS Desktop 10 and ArcPy?


texts: Script to convert line segments to equally spaced points based the duration between points as specified in another field

texts:
I have a single file containing several hundred line segments created through connecting points collected detailing animal movement patterns.  Each segment has a "seconds" field and a "uniqueID" field.

I am trying to split each line segment (uniqueID) into equally spaced points based on the seconds of time it took the segment took to complete.

For example:  If it took said animal 60 seconds to go from point A to Point B, then that line segment would be converted to 60 equidistant points.  Each segment would have a different number of points based on that "seconds" field.


texts: Exporting from ArcGIS 10.2 to PDF for editing in Inkscape-Where are the layers?

texts:
I'm trying to export a simple street map with markers from ArcGIS 10.2 to PDF for editing in Inkspace (Note: I've tried this with exports to .ai and SVG as well with results similar to those recounted below). I'm somewhat new to Inkscape having mostly used AI for polishing, but would like to give this FOSS option a try.

The odd thing I'm coming across is that while the PDF loads well in Inkscape, the layers tab either shows no layers or one for the whole map. Either way there is no listing of the different items (streets, labels, markers, etc.) in the map, nor groups. I can select items in the map and manipulate them individually, but as far as I can tell, that's about it. It seems a lot of people are using Inkscape so I'm guessing I'm just overlooking something here, but how do I access the individual elements in the map from that layers tab?

Thanks!


texts: How to get the file location of rasters registered with the PostGIS extension?

texts:
Using raster2pgsql I am able to load a NetCDF file into PostgreSQL using the PostGIS 3.1 extension. And looking at the documentation I can see that it's possible to register a raster file instead of copying the data to PostGIS using the --register flag (https://postgis.net/docs/using_raster_dataman.html):
-R, --register
Register the raster as a filesystem (out-db) raster.

Only the metadata of the raster and path location to the raster is stored in the database (not the pixels).

Then I see that there is a 40KB table in PostgreSQL that corresponds to the raster file. How can I get the file location of this raster using SQL? Surely since Postgres links to the file I can find out where the file is stored on the disk?


texts: Removing elevation trend over sloped surfaces?

texts:
Is there a known methodology for removing the trend from a DEM surface?

For example, if I had a DEM of a river, that river would flow from high elevation (upstream) to lower elevations (downstream). Even if the distance from the top of the water surface to the bottom of the river bed were the same across the whole river, the downstream portion of the raster would still have lower elevation values in the DEM. Is it possible to even out these values so that areas of equal depth also possess equal elevations?

When removing this trend it would be ideal to preserve the overall characteristics of the DEM, such as dips in the riverbed.


texts: What is the max. length of string for GP output?

texts: I have an XML as the output of Geo-processing tool. I am keeping it as a file right now but another developer wants it as return value. So what is the maximum possible length of this String output? Is this variable type or fixed type based on MS terminology?- Thanks


texts: Kaspersky Dangerous Object when installing QGIS 2.18

texts:
When I install QGIS 2.14 and 2.18, Kaspersky gives me a Dangerous Object message.  Below is a screen shot of the quarantine.

Is this something I need to worry about?


texts: Is there a way to save images with none stretching in ArcGIS

texts: I have some satellite images, when I imported these in ArcGIS and went to properties and then symbology tab, I saw stretch type was set to standard deviation. I want to save these images with no stretch as I import these images into ArcMap, to display the images without stretching , how can I do so?


texts: What are the different tools / libraries available for 3D map generation?

texts:
Like in OpenLayers 3- Cesium integration, OpenLayers3 handles the features and functionalities pertaining to the map while the Cesium counterpart handles the 3D view or the globe view.
Likewise, what are the other tools or libraries available for the 3D map generation , standalone libraries (like Cesiumjs) or integrated packages (like OL3-Cesium). I have to do kind of a functionality comparison among the commercial tools. Can anyone give me some names? It'll be even better if you can suggest a few of both commercial and non commercial.


texts: Mean coordinate issue with multipart polygons

texts:
I've been using Mean Coordinates in QGIS to find weighted mean centers based on different attributes for the counties in a singular state of the US. However I have come across a problem in which the calculation is affected by the number of polygons within a shape.

For example in the attached image I have been trying to simply find the geographical mean center of the state of Maine (no weights). The center should be further north but because the islands towards the south act as attraction factors the star, which is the result, finds its way towards the south.

Is there anyway to fix/work-around this issue? Do I just have to delete all islands?


texts:
Projection is in NAD 83- units of measure in meters.

I have major roads NC vector data overlaid on Forsyth county in North Carolina.

I would like to create a cost surface distance map using this major roads field. I converted the major roads vector to raster format using measured distance of each road in miles as the field value.

Next, I converted the major roads raster from miles to meters using calculator (divided by 0.0006213). To find the cost surface, I input Forsyth county as source data, major roads raster as input cost raster, and I received an output cost map with unrealistically high values. If I can get help to find the correct cost map that would be great.


texts: Clipping two layers to get contours underneath one layer

texts: I have a contour layer and another of the area I want to trim the contours down too.  I have them in the same CRS and i've tried intersection and clip under geoprocessing.  I can't seem to get any output from these tools. The clue is the boundary I want and the purple is what I want trimmed.


texts: Error in MMQ2GIS what is going wrong here?

texts:
To make sure things go efficiently in this green service community I created a skeleton shapefile or two, based on copy's they could give me. I created structure/style's settings by creating new attributes and value map's. I also save this locally so I can load the style in the new shapefile.

Because they stil need to get working, there shapefile contains new objects but do not have the same structure/shapesettings.

Because they are the same in spatial relation I tried the MMQGIS, to merge these two layers into one. However when I try do this an "Error" occures.

Keep in mind this in dutch: I cannot sadly translate this


texts:
I have the problem with the enabling extension SFCGAL.

CREATE EXTENSION postgis_sfcgal;

doesn't work for me even I have the PostGIS v2.2.

When I do that, it throws an error

  ERROR: could not open extension control file
  /Applications/Postgres.app/Contents/Versions/9.5/share/postgresql/extension/postgis_sfcgal.control:
  No such file or directory.


texts: Exporting georeferenced text information from a vector layer in QGIS

texts:
Translated (Google)

  "Hello! how do I export and create a new georeferenced text
  information file from a vector file that contains this information in
  Qgis 3.0.2?

  In Arcgis it is possible to export this text information and from
  there to create a new layer only with these georeferenced texts."

Title:Exportação de informações de texto georreferenciados de uma camada vetorial no Qgis

Olá! como faço para exportar e criar um novo arquivo de informações de texto georreferenciado a partir de uma arquivo vetorial que contenha essas informações no Qgis 3.0.2?

No Arcgis isso é possível exportar essas informações de texto e a partir disso criar uma nova camada somente com esses textos georreferenciados.


texts: How to obtain mean, maximum, and mininum of all points located within polygons using spatial joint in python?

texts:
I am trying to estimate summary statistics of points located inside a polygon. Each points in the point layer is associated with N attributes. My target is to summarize (say mean, min,..) the attributes of points located within a polygon and populate the attribute fields of the corresponding polygon.

I am looking for a solution using GeoPandas or other Python libraries.

import geopandas as gpd

gdf_points = gpd.read_file('/path_to_points.json')
gdf_polygon = gpd.read_file('/path_to_polygons.json')

dfsjoin = gpd.sjoin(gdf_polygon ,gdf_points)

Now, how can I summarize the stats for each attribute in the point layer and add it to the polygon shapefile? Which function can I use?

I am looking for something that is functionally equivalent to ESRI ArcGIS SpatialJoin_analysis with fieldmappings


texts: Performing bounding box query in PostGIS?

texts:
I have a PostgreSQL table, with almost 2 million rows, with a long-lat coordinates field in the form POINT(-73.4938 33.2405).

Supposing there's a geospatial index on that field, what's the most efficient, fastest way to select all the rows within an arbitrary bounding box?

The box is like SW long-lat: -74.0042 40.7688, NE long-lat: -73.8809 40.7984.


texts:
A CSV file contains the coordinates of all polygons which refers a particular crop in a specific time period.
Eg:

How I can load this csv into PostGIS..
I tried with CREATE TABLE
But I don't know which variable is used to represent geo column...ie, for representing spatial coordinates.


texts: Choosing standard parallels / central meridian

texts:
I have to choose a equal area projection for the country of Columbia and modify the central meridian and two standard parallels. I am using the South America Albers Equal Area conic projection.

I know how to go in a modify the parallels but how do I chose which standard parallels and central meridian to use?


texts: How to migrate from ArcGIS for Desktop to QGIS?

texts:
Our lab has about 10 projects in ArcGIS 9 (a couple are made in ArcGIS 10) that includes the following in each:

A georeferenced (WGS1984 geographical projection) raster seafloor image, representing areas of about 20m x 20m.
A layer of points across the image, each of which is hyperlinked to documents residing on the same computer.
A few layers of points (feature class) representing features in the raster image.
A few layers of polygons (feature class) representing features in the raster image.

With the above, we have been able to:

Analyse the areas outlined with the polygons and added buffers to some of them.
Answer questions like: Does this kind of fish (represented by a layer of points) always occur near a certain rock type (represented by polygons or points).
Add new seafloor images (as new layers) representing the same area, but taken at a later time, to see if the distribution of features in it have changed over time (i.e. there's 15% more fish, etc.).

In the future, we would like to add bathymetric data so that we can make a 3D view of seafloor represented by our images.

Finally the question: Can the above be successfully migrated to QGIS and HOW do I do it? (so we can save licensing fees and also use open source software, etc.)

Based on my initial look a QGIS, most of what I outlined above can be done, but can anyone confirm that?

If we do this, it needs to be done within the next couple of months, all by ourselves with no programming ability. (we don't have funds to hire someone to do the migration for us, or to add feature to QGIS).


texts: how to produce xml file in SLD to add discrete colors and color legend

texts:
I have to produce a xml file, according to Styled Layer Descriptor (SLD) (for more informations you can see: http://www.opengeospatial.org/standards/sld). The file is required by GeoServer and it has the purpose to add a color legend. How would this xml file be written, according to standard SLD ? Can you indicate me some link with examples?

Thank you for your attention.

Best regards, B.O.


texts: I am trying to add a .csv file into QGIS but the spatial reference data has not been split into lat/long, how do I get GIS to understand this sort of data, or is there anyway I can mass split up the data, as I have 10000 points to add.


texts: Compute the azimuth between two UTM32N points using Java GeoTools

texts:
I am using Java GeoTools to compute the azimuth between two UTM32N points, e.g. between POINT(496383,5419109) and POINT(528196,6070949).  Is it possible to do so without converting the points to WGS84 first?

The following code does not run because Latitude 5419109°00.0'N is out of range (±90°):

CoordinateReferenceSystem crs = CRS.decode("EPSG:32632");
GeodeticCalculator gc = new GeodeticCalculator();
gc.setStartingPosition(new DirectPosition2D(crs, p1.getX(), p1.getY()));
gc.setDestinationPosition(new DirectPosition2D(crs, p2.getX(), p2.getY()));
return gc.getAzimuth();

After a ton of trial-and-error I was able to get the code to run by converting the points to WGS84 first:

CoordinateReferenceSystem crs = CRS.decode("EPSG:32632");
MathTransform transformToWgs84 = CRS.findMathTransform(crs, DefaultGeographicCRS.WGS84);
Coordinate p1Wgs84 = JTS.transform(p1.getCoordinate(), null, transformToWgs84);
Coordinate p2Wgs84 = JTS.transform(p2.getCoordinate(), null, transformToWgs84);
GeodeticCalculator gc = new GeodeticCalculator();
gc.setStartingPosition(new DirectPosition2D(DefaultGeographicCRS.WGS84, p1Wgs84.getX(), p1Wgs84.getY()));
gc.setDestinationPosition(new DirectPosition2D(DefaultGeographicCRS.WGS84, p2Wgs84.getX(), p2Wgs84.getY()));
return gc.getAzimuth();

However, this "solution" is somewhat unsatisfactory.  Is there a way to directly compute the azimuth between two UTM32N points?


texts: Finding what intersects feature without join in ArcGIS Pro

texts:
I am new to ArcGIS Pro.

Which Geoprocessing tool do I use to find the Features in FeatureA that intersect FeatureB without prerforming a join?

I just want a selection of FeatureA - not clipped to FeatureB and not altered in any way.

Both Features are polygons.

I'm looking for a select one feature by another feature.

The Intersect Tool seems to create a new feature class.


texts: How to compute climatology of 3 days moving averages for rasters pixel by pixel?

texts:
I have one year of datasets that I want to compute the climatology of 3-day moving averages and then I subtract this climatology from the original values for each pixel.

example:

  # Set up the rasters
r1 <- r2 <- r3 <- r4 <- r5 <- r6 <- raster(nrows=10, ncols=10);
 # Populate them with some values
r1 <- setValues(r1,runif(100,min=1,max=100));
r2 <- setValues(r2,runif(100,min=1,max=100));
r3 <- setValues(r3,runif(100,min=1,max=100));
r4 <- setValues(r4,runif(100,min=1,max=100));
r5 <- setValues(r5,runif(100,min=1,max=100));
r6 <- setValues(r6,runif(100,min=1,max=100));
# Stack them
 st1 <- stack(r1,r2,r3,r4,r5,r6)

Any idea on how to do this.


texts: Curved labels with data defined coordinates in QGIS

texts:
My goal: I want to have curved labels placed in the center of a linestring. I define the label position with the expression x: $x_at( num_points( $geometry)/2) y: $y_at( num_points( $geometry)/2) to center the labels.

My problem: When I set data defined coordinates the label gets displayed horizontally even though I have set the label placement to curved (possibly due to an old QGIS bug):

My (unsatisfactory) solution: I also set an expression for the rotation (angle_at_vertex( $geometry , num_points( $geometry)/2)-90), but the results are by far not as good as having a curved label:

Does anyone have an idea how to create curved labels centered on a linestring?


texts: How can I create a Choropleth by combining a polygon GeoPanda with a point GeoPanda?

texts:
I have two GeoPanda dataframes, one with the location of shops with point coordinates in the geometry column, and one with London boroughs with polygons in the geometry column.
My usual approach with creating a choropleth is to loop through each borough and check how many points are contained in each borough, before assigning that value to the borough, for instance:
Points_Within = []
for B in borough_df:
    Points_Within.append(len(Shop_df.intersects(B.geometry)))
borough_df['Points Within'] = Points_Within

This is fine (albeit cumbersome with large dfs), however, I refuse to accept that there isn't a cleaner way of doing it.
If anyone can shed some light on a nicer way to do this I would be very appreciative!
James


texts: What is the practical use of panchromatic band in remote sensing?

texts: In Landsat and other satellite imagery, panchromatic band is also included (which normally has high spatial resolution). I want to know what is the practical use of a panchromatic image? Is it a must with other bands for LULC mapping or is it just an optional band?


texts: Adding hyphen in ArcGIS Field Calculator to merge two columns?

texts:
I am trying to build a field calculator expression to merge two columns into one. I am looking for the word that makes it a range.

Example:  I have 2 columns in elevation, From and To. I want to make them like "4100 - 4200" where I add the hyphen. How do I do that ?


texts: Specifying dependable layer names in QGIS Processing Modeller

texts:
I am creating a model in QGIS (3.10.3) Processing Modeller (Win 10[64]).
Final Objective
To create Voronoi polygons from centroids that are clipped to a higher, polygon, geometry.
Hurdle
The Voronoi clip produces slivers that either do (Voronoi-with-centroid) or do not ('merger-candidates') contain an centroid.
Chosen method to get over the hurdle
I create a UUID for every Voronoi (whether it has a centroid or not). Then identify the neighbours, and the shared boundary lengths as per:
https://spatialthoughts.com/2020/04/08/calculating-shared-border-lengths/. Then dissolve the merger-candidates into the neighbouring Voronoi-with-centroid with which it shares the longest boundary.
I am happy with the logic of this method and am confident it will produce the result I want.
Issue
However, my model (see GoogleDrive link below) is doing something odd.
The command create UUID does what it should, but additionally, adds a UUID to the layer name.
So, for the next instruction - 'Identify Neighbours' - when I enter the Expression dialog, QGIS offers all the '..._OUTPUT' from the previously completed instructions. But Create_UUID_OUTPUT doesn't appear to be the right thing to include in the "spatialthoughts" code linked to above.
My adaptation of that code is, so far:
array_to_string(array_remove_all(aggregate(
  layer:=  **## THIS IS WHAT DOESN'T WORK YET ##**,
  aggregate:='array_agg',
  expression:=$UUID,
  filter:=intersects($geometry, geometry(@parent))
),"UUID"))

If I double-click Create_UUID_OUTPUT on the top right box in Processing Modeller, @Create_UUID_OUTPUT appears in the code window on the left. When this is added into layer:= this then generates:
'''
"Eval Error:Cannot find layer with name or ID '' "
'''
So I looked at the "Map Layers" drop down on the lower right.
Sure enough, there is a Voronois_fixed_with_UUIDs layer, which is exactly what I would have expected, so, I thought, 'let's use that'.
However, when I double-click that, this generates  Calculated_11a2a0f5_618c_4b85_90e9_07c97cdf360b on the left.
So, the layer name itself has gained a UUID, and this changes every time I run the model.
I do not want that. I expect Voronois_fixed_with_UUIDs to be the layer name, for this not to change, and for it to appear, as is, in the code on the left when double-clicked on the right.
In QGIS desktop, the layer does appear with that name. Its properties may be seen at the bottom of this post.
So I have two questions:

How do I stop the Create_UUID instruction adding a UUID to the layer name?
What should I put in the line layer:=  ## THIS IS WHAT DOESN'T WORK YET ##,?

The QGIS model may be downloaded from my GoogleDrive:
https://drive.google.com/file/d/1CDJ_pd2alE8viwiPVH-qB70xwGfDozcC/view?usp=sharing
Layer properties in QGIS
(note, there's a UUID, but it's different - due to having run the model again):
Name:
Voronois_fixed_with_UUIDs

**Source:
MultiPolygon?crs=EPSG:27700&field=fid:long(0,0)&field=OBJECTID:long(0,0)&field=zone_code:string(14,0)&field=Shape_Length_1:double(0,0)&field=Shape_Area_1:double(0,0)&field=Shape_Length_12:double(0,0)&field=Shape_Area_12:double(0,0)&field=Shape_Length_12_13:double(0,0)&field=Shape_Area_12_13:double(0,0)&field=Shape_Length_12_13_14:double(0,0)&field=Shape_Area_12_13_14:double(0,0)&field=opcs_code:string(9,0)&field=Ward71Code:string(16,0)&field=Dist71Code:string(16,0)&field=Cnty71Code:string(16,0)&field=Pop71:double(0,0)&field=ORIG_FID:integer(0,0)&field=%24UUID:string(38,0)&uid={**f53eb83b-8da7-4f8f-a045-98a48f752bc6**}

Storage:
Memory storage

Comment:

Encoding:
UTF-8

Geometry:
Polygon (MultiPolygon)

CRS:
EPSG:27700 - OSGB 1936 / British National Grid - Projected

Extent:
446240.0011192215606570,92182.7952580302953720 : 453679.9979932568967342,96610.0049105565994978

Unit:
meters

Feature count:
41


texts: I intend to upload coordinates as a layer on top of the XYZ OpenStreetMap in QGIS. I cannot seem to grasp how to make them show up?

texts:
This is the code I have so far:
w = shapefile.Writer('shapefiles/test/testfile')
w.poly([[59.4171, 15.6113], [56.8809, 14.8257], [64.6431, 13.7139]])
w.field('FIRST_FLD','C','40')
w.field('SECOND_FLD','C','40')
w.record('First','Polygon')
w.close()

Now I would like to upload the shapefile as a layer on top of the XYZ OpenStreetMap. I have looked around but I do not know how to do so. Do I use a zip of the created folder to put the coordinates on the map?


texts: Scatterplot between raster and point shp

texts:
I have a TIFF containing soil moisture values in each cell and I have a point shapefile with only 50 points of field measurements.

How can I make a scatterplot of these 2 datasets? I know I can convert my raster into Excel files and do the comparison in Excel. Using Arcmap, how can I compare the 50 shapefile points with the underlying cells of the raster?


texts:
I  would like to display a map featuring both a marker and a track, the latter loaded from a GPX file on my local machine. I have been trying to use the leaflet-gpx plugin but since I am totally new to Leaflet, geodata, JavaScript...

I can't figure out why the track does not appear on the map, whereas the marker does. The GPX file is located in the same directory as my HTML file, and the GPX seems to be valid as the track is displayed when I upload it to an online GPX viewer.

Here is my code:

    <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/leaflet/1.3.1/leaflet.css" />

The above uses the example code found on the leaflet-gpx github site, which is why it does not use the latest version of leaflet. I would like to use the latest version in future, but first I'd like to get this simple example working (I've also tried converting the data to GeoJSON format and displaying that instead, but again no luck, so I guess I'm doing something fundamentally wrong).

UPDATE: excerpt from GPX file added

<?xml version="1.0" encoding="UTF-8" standalone="no" ?>
<gpx xmlns="http://www.topografix.com/GPX/1/1" creator="Best-Mountain-Artists" version="1.1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd">

  <metadata>
    <link href="http://best-mountain-artists.de">
      <text>Best-Mountain-Artists</text>
    </link>
    <time>2017-11-24T09:00:27Z</time>
    <bounds maxlat="50.935044921934605" maxlon="8.247442189604044" minlat="50.89056090451777" minlon="8.172424659132958"/>
  </metadata>   <trk>
        <name>Wildkatzen-Walderlebnis Winterstein </name>
        <trkseg>
            <trkpt lat="50.341790" lon="8.673659">
                <ele>337.19223</ele>
            </trkpt>
            <trkpt lat="50.340613" lon="8.674818">
                <ele>346.23657</ele>
            </trkpt>
            <trkpt lat="50.340092" lon="8.674088">
                <ele>357.13358</ele>
            </trkpt>
            <trkpt lat="50.340448" lon="8.667651">
                <ele>361.92911</ele>
            </trkpt>
...
            <trkpt lat="50.341755" lon="8.673716">
                <ele>337.40000</ele>
            </trkpt>
        </trkseg>
    </trk>
</gpx>


texts:
I am trying to read a KML file in R using this previous guide. The file is available here as a KMZ; following suggestions from other help files, I converted it to a KML by opening it in Google Earth, right-clicking the data layer, clicking 'save as', and saving as a .kml.

Here is what I've done so far:

import <- ogrListLayers("Burrows/Burrows_et_al_Nature_traj_ocean_NH1.kml")
# [1] "SST"       "SST_start" "SST_1"
attr(import, "driver")
# [1] "LIBKML"
attr(import, "nlayers")
# [1] 3
foo = readOGR("Burrows/Burrows_et_al_Nature_traj_ocean_NH1.kml", "SST")
# no features foundError in readOGR("Burrows/Burrows_et_al_Nature_traj_ocean_NH1.kml", "SST") :
# no features found

Because the output of ogrListLayers() makes it clear that the path name is correct and R can access the .kml file well enough to list the layer names, I'm confused why readOGR() isn't working. What am I missing?


texts: Using Erase tool for replacement ArcGIS for Desktop?

texts:
I have two layers (polygon type) and i want to erase the red polygon from the pale blue polygon, so i will get another polygon just without the red polygon feature.
I'm work with desktop license, so i can't use the "Erase" tool. I need ArcGIS solution. (I can do it in QGIS).


texts: Addressing just the outermost cells in a raster?

texts:
I have a DEM cropped to the extends of a shapefile with different glaciers (see here for data: https://drive.google.com/open?id=1ERFdsqDGLH1a_FbxwawE_gPm0Au0Q9vT).

DEM <- raster("***Glacier_Clip1.tif")
glaciers <- readOGR("***/01_Clipped_Data", "Glacier_Clip_Polygon")

glaciers_sf <- st_as_sf(glaciers)
glaciers_comb <- st_union(glaciers_sf)
glaciers_comb <- st_cast(glaciers_comb, "POLYGON")
glaciers_large <- glaciers_comb[st_area(glaciers_comb) > units::set_units(50000,"m^2")]
glaciers_large_sp <- as(glaciers_large, Class =  "Spatial")
glaciers_ras <- mask(DEM,glaciers_large_sp)

For further calculations I need to differentiate between all the "inner cells" and the "inner boundary cells" inside the glacier directly at the border and the "outer boundary cells" directly outside (those are not part of the cropped raster, so I somehow need to get them from the original DEM).

This image shows the three cell types (light blue, grey and yellow) I need to address somehow. Is it possible to get them via extract? Can anyone help me please?


texts: FME error: rror connecting to instance 'sde:sqlserver:Server' for database 'database' on server '!not_used!'

texts:
We are facing the error  SDE_ERROR_CODE(-92) Bad connection handle while running the FME workspace via workbench.

Full error message:

  Error connecting to instance 'sde:sqlserver:Server' for database
  'database' on server '!not_used!' for user ''.:SDE_ERROR_CODE(-493)
  Connection was attempted with an older version of SQL Server client
  communications software that is not supported.. Error connecting to
  instance 'sde:sqlserver:Server' for database 'database' on server
  '!not_used!' for user ''.:SDE_ERROR_CODE(-92) Bad connection handle.

Note: SQL native Client is SQL server 2012 native client and Microsoft system CLR types for SQL Server 2012.

We are facing this after we upgraded to FME 2016.1.2.1.
The workspace can run successfully in 2014 FME version.

Kindly suggest how to resolve this.


texts: QGIS Helmert transformation plug-in: "Missing metadata file" error

texts: I want to convert the integer coordinates with the x and y coordinates to the coordinates of ED50. When I install the QGIS 2dhelmert plug-in, I get a 'missing metadata file' error. How do I solve this problem? Is there another method?


texts: Avoiding dbo.arcsde LOCK issue while exporting to ArcSDE?

texts:
I had a spatial db connection open in arcmap. While exporting to SDE with a console app I get an error in the console:

  Lock request conflicts with an established lock [DBO].

Does this just mean someone is using the same ArcSDE connection or is it the feature dataset?

Of course when I closed arcmap, and fired the console app again, my feature class exported to ArcSDE successfully.

Can I kill any locks programmatically if I am DBO (when other users are using the same FD)?

If that's not smart, how can I avoid this lock?

How does the GIS community deal with this issue, wait till all users are off and run it at 2am?

The console app will become a batch process updating FCs on a regular basis.


texts: Is it possible to save multi ring polygons to the ArcGIS server REST API?

texts: I'm trying to insert new polygons that the user has drawn through the addFeatures REST endpoint and I'm getting an error when the user draws a bow-tie and tries to insert it.  Is this simply not supported or am I missing some geoprocessing step?


texts: How to overlay vector grid to keep track of areas I have reviewed:ArcMap

texts: I need to review large areas at the street level (in ArcMap) to make sure certain features are up to date. Is there a way to overlay a grid so I can "check off" areas I have already covered? I need the grid to show up in data view and layout view.


texts: I am doing some research on consumer grade GPS accuracy and precision. I took 5 readings for 20 known locations (100 readings). Am I correct in assuming that the standard deviation between the 5 readings at each location would indicate precision, whereas the average distance from the five points to its known location would indicate accuracy?


texts: Importing geometry field without clicking run in google earth engine

texts:
Hi I create a google earth engine code. I want to run this code without clicking run. So I create a button that runs my all function. Also I get my date inputs with Textbox.
However; the issue is when I select the geometry field and click my button that runs my functions, the code says you must select geometry. To solve this I must click to run button after the geometry selection and then the code works well. Is there an any method for importing geometry without clicking run or can I create a button that importing selecting geometry into code?


texts:
I am cross-checking an attribute table and clicking on features one-by-one. Is there a tool or key to jump to the next feature (by feature ID), instead of having to click in the gray cell to the left of a record, like a "next" shortcut?

I am using ArcGIS 10.3.1.


texts: How can I make a map-marker's popupContent be a random phrase?

texts:
I'm using leaflet.js to make a map and I need help with the javascript required to make a map-marker generate a random phrase with each click.

I've used this code to generate random text:

function randomText() {
var items = new Array(
            "Ow",
            "Stop poking me",
            "Stop",
            "Leave me alone",
            "Don't you have anything better to do?",
            "You're hurting me",
            "Go away"
);

var i = Math.floor(Math.random() * items.length);

return items[i];}

But when I try to put it in popup content it's blank.

I tried

.bindPopup(document.write(randomText()));

My current gist has .bindPopup(randomText()); because it is, at least, not blank and  generates a random phrase with each page load.

Here is my gist https://gist.github.com/joshuahouston/800cd59aab26c264f3bc

And the working map with where I'm currently stuck http://bl.ocks.org/joshuahouston/800cd59aab26c264f3bc


texts:
I am using QGIS 2.8 and am trying to create/load an SCP file in which to classify a Sentinel 2 image. It won't let me create an SCP file and keeps asking if I want to add missing fields to the shape file. If I answer yes or no, it still won't create an SCP file. I've tried uninstalling and reinstalling the plug-in but it doesn't help.

I tried to update to to 3.4 but that didn't work at all so I went back to 2.8. The SCP file is the training file to put the ROIs into. Without it I can't create ROIs.


texts: Have Atlas sorting and filtering controls been removed in 2.4?

texts:
In this diagram

the atlas generation tab of the print composer has several controls that I don't see in my version (2.4) of QGIS.  In particular the feature sorting and feature filtering, both of which sound like they might solve a problem I am having with Atlas generation see :How do I control rendering order of items in composer?.

Is this a screen shot from an older version and these controls have been removed?


texts: underwater mosaic with overlapping photos

texts:
What is the best approach to creating a mosaic with underwater photos.

The photos are black and white. They suffer from non-uniform illumination, i.e. the middle of the photo is bright and the edges are black. My attempts creating a mosaic with such photos result in a mosaic with stripes marks.

An example of an ideal result is attached.

Left : before and right : after

AUV, underwater, illumination-correction


texts: Openlayers WMTS GET request with Basic Authentication returns 401 Error

texts:
We are trying to publish WMTS with basic authentication from Geoserver.

Our addwmts() function works without a problem when there is no basic authentication.

We have stored our tiles in Geoserver. We have placed a Rest API between Geoserver and the Web interface.

We are trying to publish WMTS tiles that we get from Rest API with basic authentication onto the web interface.

Our Rest API returns us a .xml document and we want to read this .xml to visualize it on our web interface.

Relevant parts of our JS and html codes are below. We have a button at our index.html file that is triggering addwmts() function.

When we run below code we get

We can not send authorization header for our tiles. How can we fix that 401 problem ?

We have checked most of the links including the following :

OpenLayers 2.12 and http basic authentication woes

Openlayers 3 getting 401 error from geoserver with secured layers

<!DOCTYPE html>

    <meta charset="utf-8">
    <meta name="viewport" content="initial-scale=1.0, user-scalable=no, width=device-width">
    <title>trial</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v4.6.5/css/ol.css" type="text/css">
    <link href="asdf.css" rel="stylesheet">

    <button onclick="addwmts()">WMTS</button>

  var parser = new ol.format.WMTSCapabilities();

  function addwmts() {

var username = 'myUserName';
var password = 'myPass';

function make_base_auth(user, password) {
    var tok = user + ':' + password;
    var hash = btoa(tok);
    return 'Basic ' + hash;
}
$.ajax({
    type: "GET",
    url: "http://xxx.xxx.xxx.xx./xxxxxx/xxx/getCapabilities",
    dataType: 'text',
    async: true,
    data: '<>',
    crossDomain:'true',
    beforeSend: function(xhr) {
        xhr.setRequestHeader('Authorization', make_base_auth(username, password));
    },

    success: function(text) {
        let result = parser.read(text);
        let options1 = ol.source.WMTS.optionsFromCapabilities(result, {
            layer: 'myIdentifier',
            matrixSet: 'EPSG:4326',
            format: 'image/png',
            projection: 'EPSG:4326',
            crossOrigin: 'use-credentials'
        });
        var overlay1 = new ol.layer.Tile({
            title: 'myLayerTitle',
            preload: 4,
            opacity: 1,
            visible: true,
            source: new ol.source.WMTS((options1))
        });
        map.addLayer(overlay1);
    }
}); };


texts: How to use Twitter Typeahead with GeoServer WFS?

texts:
I'm trying to configure the Twitter Typeahead library to work with a GeoServer WFS layer, so I can create an auto-complete feature search function.

Typeahead has the Bloodhound suggestion engine, which can work with a remote data source. There's a good working example at http://jsfiddle.net/Fresh/bbzt9hcr/ where the transform function is used to convert the nested JSON response into an array which is displayed in the Typeahead dropdown.

I believe that it's necessary to use JSONP and a callback function in order to avoid cross-origin issues in GeoServer (please correct me if I'm wrong) so I'm using the following code to achieve this:

    var parameters = {
        format_options: "callback:typeaheadCallback",
        service: "WFS",
        version: "1.0.0",
        request: "GetFeature",
        typeName: "repairernames",
        outputFormat: "text/javascript",
        propertyName: "repairername",
        CQL_FILTER: "repairername like 'QUERY%'"
    }
    var repairers = new Bloodhound({
      datumTokenizer: Bloodhound.tokenizers.obj.whitespace('value'),
      queryTokenizer: Bloodhound.tokenizers.whitespace,
      remote: {
        wildcard: '%QUERY',
        url: 'http://localhost/geoserver/blah/wfs?' + jQuery.param(parameters),
        transform: function (result) {
            return $.map(result.features, function (feature) {
              return {
                value: feature.properties.repairername
              };
          });
        },
        prepare: function(query, settings) {
          settings.dataType = "jsonp"
          settings.url = settings.url.replace('QUERY', query)
          return settings;
        }
      }
    });
    function typeaheadCallback(json){
        console.log("Typeahead callback response: ", json);
        return json;
    }

The change compared to the working demo is the new prepare function (described here) which tells Bloodhound/Typeahead to use JSONP, and the callback typeaheadCallback which is specified in the parameters.

Via the debugger I can see that the request to GeoServer is successfully implemented, and the console.log in the typeaheadCallback shows that valid values are returned as JSON.

The problem is that the transform function is never called, so the returned JSON object from GeoServer isn't transformed, and therefore isn't sent to the Typeahead dropdown menu, which therefore doesn't display results.

How can I configure Typeahead and Bloodhound to work with a GeoServer WFS layer?


texts: Convert 5-digit string number into a readable date in CartoDB

texts:
Based on a couple of assumptions here

1) I have about 6 clean formatted date columns imported from Excel that CartoDB converted into 5-digit strings. I cannot select 'date' in the column label.

2) Assumption two is that it's realistic and logical to change the dates so they're readable. (that is, convert them via SQL to eg YYDDMM).

I got help from HQ in the form of to convert the columns use:

SELECT to_date(date_column::text, 'MDDYY') AS fixed_date

Since I don't know any other way to run SQL commands against my db, I'm trying to use the SQL API. After trying everything I could and from Google I just don't see why the above query would be useful. That is, was it perhaps offered as a per-map command I should put on relevant layers? So my approach is to try to run from REST but, shouldn't that then be a POST or UPDATE? Here's my latest:

https://myusername.cartodb.com/api/v2/sql?q=SELECT FROM myTableName to_date(date_column::text, 'MDDYY') AS fixed_date&myApiKey

At the end of the day I just want to change many date columns to readable dates that people can filter from and that anyone can read.


texts: Updating 'sum line lengths' polygon layer when line layer is updated

texts: After using sum line lengths, is there an automated way to update the resulting polygon layer when more lines are added without needing to delete the layer and re-run the sum line lengths algorithm?


texts: QGIS - Semi-Automatic Classification Plugin v.4.9.5 - When trying to create ROI getting Error[43]: Attribute table error

texts:
when I try to create a ROI within the Plugin: Semi-Automatic Classification Plugin, I am getting the Error: Attribute table error. Which means it is not creating MC ID / MC Info or C ID / C Info in the ROI list. So no attributes are created, but the weird thing is the signature list gets the correct informations.

So far I tried reeinstalling it, which didn't work. Any hints how to solve that problem?

Hope the picture makes my point more clear:

Still trying to solve this problem. I don't understand why it's not copying any attributes to my created ROIs in the ROI list.

When I try to delete the created ROI, following error appears:
AttributeError: 'NoneType' object has no attribute 'text'


texts:
I have java-script menu that is getting behind openlayers OSM layer (only after it loads).
I've tried giving z-index:9999 to the menu div but it didn't work.
Any idea ?


texts: I have my own map matching algorithm and I want to compare it's accuracy to some other algorithms. Is there a free trace dataset that has a confirmed match for each trace (ground truth), so I can perform the comparison? Even more, I would be interested in a whole benchmark/test suite for this purpose.


texts:
I'm trying to create a new layer, in which the buffers are displayed. I tried this code, but it does not work. Where is the error?

layer = iface.activeLayer()
buff_layer = QgsVectorLayer("Polygon?crs=epsg:4326","buffer", "memory")
bufferDist = 10
for feat in layer.getFeatures():
  geom = feat.geometry()
  buff = geom.buffer(bufferDist,5)
  feat.setGeometry(buff)
  buff_layer.addFeature(feat)
QgsProject.instance().addMapLayers([buff_layer])


texts: Choropleth map from voronoi map from a set of coordinates (d3+leaflet)

texts:
I have a csv file with a set of coordinates with the following structure:
id,name,latitude,longitude
(24 values).

Those 24 points represent coordinates in a city. I would like to get a voronoi map with those coordinates as centroids, and then colour the polygons as if it were a choropleth map. The result of this will be a map to be published on a web dashboard, so I will be using js-only frameworks.

I've been following this tutorial and this example (leaflet+voronoi), but it's hard to find:

Tutorials with current/modern versions of d3 and/or leaflet.
Tutorials or articles on how to do the mapping process I need (coordinates>voronoi>choropleth).

Currently I am following the d3.voronoi section in the d3 API docs in order to understand how to generate the regions/polygons, but I don't know how to process those in order to generate the choropleth map using leaflet and mapbox.


texts:
I have for now just used the standard Get Feature Template in the documentation. The first example should make all the returned data format into bullets instead of little tables. For some reason my middle returned feature still shows as a table. And it appears Geoserver is not using the template for this layer, visible from the 333 I typed at the end of the feature.type. Any suggestions to why GeoServer is not formatting the middle table? It must be a setting I am guessing.

All the info returned is coming from PostGIS views, the GeoServer has the same access rights to each view.

Update
I had a play, and if I remove the content file from geoserver_data/templates and place it in geoserver_data/workspaces/<workpacename>/<store>/<layer> the template loads again for the first layer. However if I do the same for map_abundace_overview again the template is not applied. I did a search for tfl files and found the one i created and a couple under medford, which i renamed. Any suggestions?

<#list features as feature>

Type333: ${type.name} (id: ${feature.fid}):

  <#list feature.attributes as attribute>
  <#if !attribute.isGeometry>

${attribute.name}: ${attribute.value}

  </#if>
  </#list>

</#list>


texts: GDAL Python - How I can determine the PROJ version used

texts:
As PROJ recent major update introduces importants change in the API, I need in my script to check what is the PROJ version compiled with the current GDAL installation of the user. I know I can get the GDAL version with osgeo.gdal.VersionInfo() but it does not tell me what is the PROJ version.

Any suggestions?


texts:
I have a problem with PostGIS and a connection to QGIS. I successfully connected to a database on a server but when I try to add a PostGIS Table, I am getting warning symbols next to all my tables. I have permission to all the tables and
I am using the latest version of QGIS.

Any hints?


texts: How to delete MDB file after using UpdateCursor?

texts:
I am updating MDB file with update cursor after updating complete I want change the file extension MDB to other. It shows it's in use by another process. Find below code.

arcpy.env.workspace=ConvertedMdbFile
Styletablelist=arcpy.ListTables("","")
for table in Styletablelist:
    if str(table)=='Fill Symbols' or  str(table)=='Line Symbols' or str(table)=='Marker Symbols':
       with arcpy.da.UpdateCursor(table,["Name","Category"]) as Tablecursor:
          for row in Tablecursor:
             namelist=row[0].split(".sde.")
             if str(row[1])=="Layer" and len(namelist)==1:
                 styleName=row[0]

                 row[0]=prefixname+styleName
                 Tablecursor.updateRow(row)
del Tablecursor
del table
del Styletablelist
OriginalMDBFile=str(os.path.splitext(INPUTMDBFILE)[0])
os.rename(ConvertedMdbFile , OriginalMDBFile+ '.style')

In above os.rename through error

  [Error 32] The process cannot access the file because it is being used by another process

How can I resolve this?


texts: Speeding up GeoServer response - splitting the main PostgreSQL table into multiple sub tables instead of using views would help?

texts:
I am sorry if the title might sound silly, but I will try to explain the problem and why that title seemed the right to me.
I have an OpenLayers application which display different WMS Image Layers and make query against the relative WFS services coming from GeoServer.
These layers are published from different PostgreSQL views. All the views are generated against a join with a main attribute table and a main geographic table through common feature ids.
The main attribute table is getting big (right now it is 2.5 millions rows), and the response of the application in serving the WMS and other operations like using cql_filter is getting slower.
I cannot create a tiled WMS as my data change frequently and data is served with a filter on dates, so that I would have to create tiles for each single day I have in the views.
Also, I took care of serving the layers with the same coordinate system as my data not to slow down the process due to reprojection.
I tried to create SQL Views in GeoServer to further decrease the size of each view in order to see if this would have improved performance, but I still have slow responses, and I guess (I might be wrong though), that this is due to the data eventually coming from the same big main table under the hood.
So, my idea was to split the main big attribute table into smaller tables, and create my views against these.
Would this be a valuable path or do you have better suggestions?
I can provide more information if needed.


texts:
I have a large dataset in a PostGIS database. I am viewing it via a web interface that loads the data as GeoJSON. There is too much data to give the client all at once (10s of MBs worth...) so I would like to get an area of data within the bounding box formed by the edges of the window (the dotted line in the image below). Getting the coordinates for points f and g is easy.

Question 1: Is this a good way of doing things? Should I be thinking of caching them as tiles or is this method likely to be efficient enough?

Question 2: How do I retrieve just the data within this bounding box?

Question 3: If a shape overlaps the edge of the bounding box (e.g. shape A below) is there a simple way to crop it as it is queried like in the second image?


texts: Software for matching spectral end members with spectral library

texts: I am working with Moon Mineralogy mapper(M3) - hyper-spectral data and i am trying to match the endmembers obtained from M3 with RELAB spectral library. I have tried spectral analyst tool in ENVI but failed to match spectra due to spectral mixing of pixels. I am looking for a software or module or technique that matches each endmember with various possible materials from spectral library.


texts: Getting full path of empty sde feature datasets (Children property from arcpy.Describe object won't work))

texts:
My goal is to copy all the feature classes in a file geodatabase into an sde geodatabase, skipping those elements which already exist in the output sde gdb. Both geodatabases have the same feature dataset structure and the sde geodatabase may have empty feature datasets.

To copy feature classes into feature datasets I need the full path of the output feature dataset where I want to copy each. For example, for UrbanData feature dataset it would be sqldb.DBO.UrbanData.
To get that information I use this piece of code:

#input parameters
from_bd = r"C:\GISData\Buildings.gdb" #file gdb
to_bd =r"Database Connections\city_dat.sde" #sde connection file

arcpy.env.workspace = from_bd
sdebd = to_bd
datasetList = arcpy.ListDatasets()
desc = arcpy.Describe(sdebd)

for dataset in datasetList:
    for child in desc.children:
        if child.dataType == "FeatureDataset" and child.name == dataset:
            fdset =child.name #the feature dataset full path

And I use it to copy each feature class to the output sde gdb:

    for fc in arcpy.ListFeatureClasses('','',dataset):
        if arcpy.Exists(sdebd + os.sep + fc)==False:
            arcpy.Copy_management(fc, sdebd + os.sep + fdset + os.sep + new_fc)

To my surprise the Children property from arcpy.Describe object won't return empty feature datasets.
There's no mention of this behaviour in the Describe Object documentation.
So I need a workaround and I'm a bit stuck. Any ideas?


texts: Counting number of fields in table (or feature class) using ArcPy?

texts:
ArcPy contains a ListFields method, which:

  Lists the fields in a feature class, shapefile, or table in a
  specified dataset.

The resulting object has a count method, but I can't see how to return the field count from this (shouldn't this return an integer?).

Aside from iterating through the fields and incrementing a counter, how can I return the field count?

import arcpy
fields = arcpy.ListFields(<table>)
count = fields.count

[Dbg]>>> count
<built-in method count of list object at 0x0EDEABE8>


texts: How to Aggregate and dissolve MySQL geometry

texts:
On CentOS, I have a MySQL 5 database with a MyISAM table that has a geometry column with a spatial index, containing US county boundaries.

Trying to combine/aggregate/dissolve some of these into as few polygons as possible, to create minimum KML/KMZ files for public use, using the ogr2ogr just built with the GDAL 1.10.1 sources.

INFO: Open of `MySQL:mydb,user=root,password=password'
      using driver `MySQL' successful.

Layer name: counties
Geometry: Unknown (any)
Feature Count: 3234
Extent: (-178.309480, -14.601813) - (145.909365, 71.441059)
Layer SRS WKT:
(unknown)
Geometry Column = geo
county_fips: Integer (0.0)
state_fips: Integer (0.0)
state_abrv: String (100.0)
county_name: String (255.0)

@afalciano provided the syntax for combining the data for ESRI Shapefiles, but I have virtually no GDAL/ogr experience, and am having problems translating this to work with MySQL.

Ideally, I'd output as KML, but have had issues with that, so would be happy to get the output in WKT or anything Perl could be made to parse, or into another MySQL table with geometry column.

Can you provide any tips on ogr2ogr tutorials (other than the online driver doc pages)?

I need to do this programmatically, on demand, on a headless CentOS web server.


texts: Georeferencing window won't open QGIS 2.4

texts:
I had the georeferencer window open in QGIS 2.4 while project was running, I accidentally closed the Qgs document with the georeferencing window open (worked absolutely fine till then). When I opened QGIS again and clicked the georeferencing plugin, the window doesn't open.

I tried closing and opening the Qgs doc, deleting temp files, opening QGIS 2.2, restarting computer, but the problem persists. I don't see any error logs in QGIS when I try to open the window.

Any help or possible error checklist is welcomed.

Update
I've tried a quick re-install of QGIS 2.4, still no luck. looks like a complete purge is needed.

Update 2
I created a new user on my computer, it already has QGIS installed, and the georeferencer plugin works fine.


texts: Importing XMP Metadata into attributes fields in QGIS

texts:
I would like to import EXIF & metadatas of JPG images into QGIS.
For each image, a vector were created by 'ImportPhotos' pluggin, it adds some EXIF into attributes fields but some importants EXIF informations are missing.

In fact, I need 3 more fields. The XMP tags  I need are: GimbalYawDegree, FlightYawDegree, RelativeAltitude value are real. That's why I've created 3 new columns (attributes fields) named: GimbalYawDegree, FlightYawDegree, RelativeAltitude.

I would like to create a python script to run a batch for every feature (img imported by ImportPhotos) of a layer. The task would :

import the XMP datas of images stored at the path located into the attribute field named "Path"

import PIL.ExifTags
exif = {
    PIL.ExifTags.TAGS[k]: v
    for k, v in img._getexif().items()
    if k in PIL.ExifTags.TAGS
}

This code should read the EXIF but this is my first experience in Python and I'm not able to finish the job.

there few samples of images including GimbalYawDegree, FlightYawDegree, RelativeAltitude tags
https://we.tl/t-BjPuzBGMRc


texts: Show labels for features that fall outside of inset data frame

texts:
I'm working with a scripted Mapbook in ArcMap; There are several pages that require an inset map to show detail of a cluster of points in the main data frame.  The cluster of points comes from one layer that has labels turned on.  Only labels for features that are not in the inset data-frame should be turned on.

In the example below, only the label for SD01 would show in the lower-right, as it does not fall within the inset data frame.

My initial thoughts are to script a select-by-location (selecting points within the inset data-frame extent), building a query based off of the selection results, and applying that query as a label definition query in the main data frame.  This seems overly complicated for the task at hand.


texts: Opening shapefile created by ArcMap in QGIS

texts:
I've recently started working with QGIS after my ArcMap student licence expired. I'm trying to open files in a geodatabase created in ArcMap but most of the files in it are missing (.shp). The program simply doesn't recognize them (though a browse through the file itself shows they are there - not accidentally deleted).
If this is a compatibility problem, then how can I get access to the files in QGIS without using ArcMap?
I now also notice that rasters are not opening and that the files are locked


texts: Building heights above DEM through field calculator

texts:
I have a series of building polygons which I would like to convert to a raster illustrating their height. I have the building heights but not their elevation above sea level datum. I have the underlying topographical data in a DEM. So, I was wondering if someone could show me an algorithm in the Field Calculator that will allow me to combine the two to arrive at a building elevation above sea level.

(Apologies if this is an easy question, I am a newbie to GIS, or has been asked before but I couldn't find it in the previous questions)

Yours in hope


texts: One-sided buffer command line doesn't work?

texts:
I'm trying to create a one-sided buffer from a polyline shapefile using ogr2ogr command line and a virtual file. However, it only creates a .dbf file and not the shapefiles.
Here's the code:

ogr2ogr -sql "SELECT ST_Buffer(A.geometry,50, 'right') AS geometry, A.* FROM Polyline A" -dialect SQLITE Buffer.shp input.vrt

I have to mention that this code works fine when not putting the 'right' or 'left' parameters.

Any ideas?


texts: Problem with the ST_Contains function in PostGIS on the comparison query


texts: Simplify polyline with the intersection point kept

texts:
I have a some line and polygon features which consist of too many points, so I tried to simplify them.

I found the "Ramer–Douglas–Peucker algorithm".

I think it will work well for a single feature (line or polygon), but in my case, the feature are connected with others like this:

Note the point inside the red circle is the intersection point of the two line. This point must not be removed during the simplify process.

So I wonder if there is any out of box solution?

BTW, we have ArcMap Desktop 10. But we prefer to make it by program rather than manually since there are almost 20+ layers.


texts: Handling Tiff+tfw with High and low value using GDAL and C#?

texts:
I am converting tif+tfw to Geotiff using GDAL C#, below is the sample code. I can able to handle the color images with Band 1, 2 and 3 tif images, but for the Black and white images which had no band values, it contains high and low values, it throws exception.

WriteEncodedTile/Strip() failed.

OSGeo.GDAL.Driver.CreateCopy(String utf8_path, Dataset src, Int32 strict, String[] options, GDALProgressFuncDelegate callback, String callback_data)

GdalConfiguration.ConfigureGdal();
                String[] options = { "COMPRESS=DEFLATE", "PREDICTOR=2" };

                string fileNamewithoutext = Path.GetFileNameWithoutExtension(fileSrcPath);

                //Gdal.AllRegister();
                OSGeo.GDAL.Driver srcDrv = Gdal.GetDriverByName("GTiff");
                outputfile = txtSelectOutput.Text + "\\" + fileNamewithoutext + "_Georefereced.tif";
                Dataset srcDs = Gdal.Open(fileSrcPath, Access.GA_ReadOnly);
                Dataset dstDs = srcDrv.CreateCopy(outputfile, srcDs,1, options, null, null);


texts:
I'm having issues adding a layer. I'm very new to the application, so most likely it's just some basic thing I'm doing wrong.
I'm trying to add a Layer (https://data.geo.admin.ch/ch.swisstopo.geologie-geotechnik-gk500-lithologie_hauptgruppen/) onto the Swiss map (https://shop.swisstopo.admin.ch/en/products/maps/digital_maps/vector/smv500).
The result I want to further work with should look like this: https://s.geo.admin.ch/8a83336eeb
I tried several Layer-Adding-Options and stuff, but more like trail-and-error-style.


texts: What are the rectangular black markers shown in the roads of Google maps?

texts:
I'm doing some image processing to quantify the connectivity of London's roads. I'm not a geographer or map-expert and I would like to know what are the black rectangular nodes/markers in this map.

source. The map appears as shown when imported into Mathematica.


texts: What projection is this and how to use it in a GIS?

texts:
I have a request to do a map in non-standard projection:

(source: http://www4.agr.gc.ca/AAFC-AAC/display-afficher.do?id=1263246127905&lang=eng)
I couldn't manage to find it my materials. As well I wasn't able to apply such projection in ArcGIS.
Do you know what projection is it? And how to apply it in ArcGIS or Manifold or any open source GIS?


texts: ESRI Javascript API QueryTask error handling

texts:
I'm struggling to find the proper way to handle errors in a QueryTask.  Particularly if the user clicks outside the feature layer.

My current QueryTask looks at a field (bare_earth) populated with either a 1 or a 0.  I thought if I used...

if (bare_earth === 0) {
    do X;
}
else if (bare_earth === 1) {
    do y;
}
else{
    error handling;
}

However, clicks outside the feature layer return the following error:

TypeError {stack: "TypeError: Cannot read property 'attributes' of un…   at c (http://js.arcgis.com/3.9/init.js:74:436)", message: "Cannot read property 'attributes' of undefined"}
"TypeError: Cannot read property 'attributes' of undefined

Makes sense, since there are no attributes for the undefined area.  That being said, how do I go about handling this?

For frame of reference, each of my three potential outcomes (1, 0, no data) alter the value of a variable (result) that is used to populate a DIV.  For example, I want clicks outside of the feature to populate the DIV with "...outside the area...".

Thanks!


texts: Grouping polygons by nearby neighbors and by area in QGIS 3

texts:
I would like to group polygons of the same layer into groups by closer neighbors until they meet a certain area. I use QGIS 3.4. That is to say, the polygons must be united fulfilling two conditions:
1. to touch
2. That the sum of that group of polygons does not exceed more than 5 m2


texts: In QGIS heatmap howto establish a fixed maximum range value for comparisons over time?

texts:
I have several data sets that represent a type of activity within a specific geographic boundary. Each heatmap I create shows the relative activity spread nicely, but I have been frustrated in my attempts to calibrate absolute levels between each of the heatmaps. I am wanting to reveal not only the shifting in relative activity spatially, but also in terms of absolute activity overall.

For example, I'd like to create two heatmaps, say day versus night, and show how the activity moves from, say, city center to suburbia. But I also want to see if the activity is half as much or twice as much during the daytime peaks as compared to the nighttime peaks.

The data I have are a series of lat/lon pairs with a value ranging from 0 to 50, one value per dot:

Here is the second set of data, same idea, except that the maximum values only reach 30, not 50. I am trying to know how to make sure that the color for maximum in the second show that the activity not only moved, but is half the amount (quantitatively) than the first data set:

I am also having some confusion about the differences bewteen heatmap plugin and the heatmap style: how do they differ? I am on QGIS 2.8 (Gentoo Linux).  The style seems to work, but then gets corrupted, and the heatmap plugin seems kaput (it worked in 2.6...)

Here are the settings I used in 2.6 that worked (but do not at all in 2.8):

Here is a dialog from the Heatmap style I tried to use. This works in 2.8 but gets corrupted if I add other layers:

If it's not obvious yet, GIS is not my specialty.


texts: Merging two GeoJSON files for Leaflet Map in R

texts:
Does anyone know how to convert 2 GeoJSON files into 1 in R? I'm trying to build a Leaflet map for the two files together.
library(geojsonio)
library(leaflet)

dataset1 <-geojson_read("https://opendata.arcgis.com/datasets/3525b0ee6e6b427f9aab5d0a1d0a1a28_0.geojson", what = "sp")
MassData <- geojson_read("https://raw.githubusercontent.com/NewtonMAGIS/GISData/master/Massachusetts%20Town%20Boundaries/MassTowns.geojson", what = "sp")
dataset2 <- MassData
dataset2 <- subset(MassData, TOWN %in% c("BROOKLINE", "CAMBRIDGE","MEDFORD", "SOMERVILLE"))

merged <- TO_GeoJson(c(dataset1, dataset2))


texts:
I am trying to plot a feature space (scatter plot) for some SPOT images in QGIS 2.8.1. I know a histogram is available in the file properties, but I want to plot eg. Band 1 vs Band 2 to see the feature space.

The SCP plug-in  has a function which only shows the scatter plot for the created ROI. What I am looking for is to show a scatter plot (feature space) of two full bands. And if possible, also show my ROI or signatures in the full scatter plot.

I am evaluating QGIS for GIS and RS introduction for our students. We consider the feature space as an important tool for the students to understand which landscape elements the different bands are good at identifying.

In ERDAS it would look something like this:


texts: QGIS Topological Editing stops working after adding Google Maps or OSM underlay

texts:
When I work on a polygonal layer, I turn on 'Topological Editing' from the snapping toolbar, which allows moving the nodes from adjacent polygons together.

Everything works fine until I add any underlay from the QuickMapServices plugin - Google maps, OSM etc.

After any such layer is added, only individual vertices are selected. So, in order to move a node I have to move individual vertices from neighbouring polygons. The 'topological editing' button is still 'pressed' - no changes in the snapping toolbar.

The problem persists even when I remove the underlay.

How can I keep Topological Editing functional with the google/osm map enabled?


texts:
I am trying to change the CRS of one shapefile that is originally in NAD27 to WGS84.
When I proceed to make the change, I still continue to see the file in UTM coordinates...
I've tried to save the file as a new one with the and change it at the "set layer CRS" option...
Do you know another alternative?


texts: How to change polygon display/symbol properties in QGIS while digitizing?

texts:
I have two questions. It would be great if anyone help me to sort this out. Thanks in advance

1) While capturing a polygon in QGIS, how to remove the red - shade - sort of thing..... ?

2)Though the symbol of the polygon has been set as No brush - Solid line .....  While using split tool and when we select the polygon..it's get selected as a whole yellow fill. So if there is background layer which i'll be using as the reference layer to split the polygon...i can't see that one.

Split tool works even without selecting the particular polygon. But i wanna know..is it possible to set the select properties like..a simple No fill - Red Solid line ?


texts: remove unwanted areas which are not inside the boundary of second raster

texts:
1 raster layer - global map ;
2 raster layer - berlin

I cropped out 1 layer using extent of 2 layer and then overlapped these two layers however I see that 1 layer still has areas apart from berlin(1 layer). What would be the best way to remove those areas and only get the area exact as that of 2 layer?

Also, I tried using  Clipper too but somehow it is not working.


texts: Validating whether PDF file Uploaded as part of Geoprocessing Service of ArcGIS for Server?

texts:
I'm developing a Python Script that will be a Geoprocessing Service in ArcGIS Server.

The Script has a Button for Users to OPTIONALLY Upload a PDF.
If the PDF is Uploaded it gets inserted into the pack.
If no Upload I want the Code to Skip appending Page.

How would I perform the Check and pass if no Upload?


texts: Creating time series from NetCDF with multiple values to multiple points on ArcGIS Desktop?

texts:
I am using ArcGIS Desktop to read and extract netcdf files which have data on longitudes, latitudes and daily temperature over several years. I have separate shapefile with postcodes and the goal is to create a temperature timeseries for each postcode.

I am using ModelBuilder:

In my model, I initially created raster file from the netcdf file (daily_meantemp_Layer1), I then select the time dimension for each day of the month using the 'Select by Dimension' tool, then for each raster created by I extract the temperature values using the 'Extract Multi Values to Points'. The problem is that my final extracted values are all the same, I have added a table to show what I am getting

 .

It appears that my select by dimension bit of the model is selecting information from the only one or two days.

In a separate model I tried saving all the individual rasters in .gdb and extracting from them using the 'abs' command, this too gives me the same value across all days.

What am I doing wrong?

I am open to doing this using script its just I wouldnt know where to start!


texts: ArcGIS REST Server Object Extension Examples

texts: I'm looking for a few good examples of a REST server object extension (SOE) showing more than the common ESRI sample. I'd like to see some RestOperation samples as well and a few different TryGet statements (i.e.. returning integers and strings). Looking at ver 10.2 using SQL.


texts: osm2po way tag resolver: how to use allow/deny for complex selection

texts:
How to combine allow / deny in osm2po config file for complex selection?

As an example, I would like to include in my finalmask (car) all [highway=service] except the ones that DO NOT have the tag [surface=asphalt].
I would like to use this asphalt restriction only for highway=service but not for other types of highways. What should I write in my osm2po config file to achieve this?


texts: Selecting Features from the results of topology checker in QGIS

texts: Is there is any way to select or extract features from topology checker plugin in qgis?


texts:
Let's assume I have a raster image consisting of 0's and several donut-shaped rings of 1's. This raster could for example be the result of classifying a satellite image with ponds and lakes in it into two classes: 1) pond shoreline, pixel value 1 and 2) other, pixel value 0. See image for impression.

My challenge is to fill each donut hole, creating full circles with pixel value 1. For this I have access to ArcGIS with all common functionality, as well as ENVI and Python with gdal bindings. See second image for impression of the desired outcome.I guess a majority filter will work to some extent, except when the donut holes are large.

Any suggestions on transforming donuts to circles are welcome. Bonus points for making it work with unclosed donuts as well.


texts:
is there any available example online, or anybody could help me to figure out how can I use the GoogleStreetViewPanel widget in the opengeo suit sdk? I tought it will work like this:

{
        xtype: "gxp_googlestreetviewpanel"
}

but this code generate me an error about missing lan/lot...
I could manage to add other widgets, like scaleoverlay.


texts:
I'm new with MapServer and I followed various tutorials to display a shapefile in MapServer. Now I want to use my shapefile but I have this error:

  loadMapInternal(): General error message. Given map extent is invalid.
  Check that it is in the form: minx, miny, maxx, maxy

my extent field is:

EXTENT 4567164.94 4570286.94 614923.216 618804.216

These information come from ArcGIS (my colleague opened the shapefile in ArcGIS).

Can someone help me to understand what is wrong?


texts:
QGIS 2.x never had an option to export as a layered pdf. Instead there was a cumbersome workaround, involving exporting as an svg, importing into Scribus, and manually adding individual features to separate layers. All that is covered under this question:  Can QGIS preserve layers when exporting to PDF?

Is there a more efficient way to create a layered pdf in QGIS 3.0?


texts: ogr2ogr DXF to PostgreSQL some LINESTRING coordinates mirrored to the negative longitude

texts:
I have been using ogr2ogr GDAL tool to migrate geospatial data from DXF file to PotgresQL database without much problems. But now there is a special problem I can not fix.

I have a DXF file with LINESTRINGs that when I import them to the PostgreSQL database some of the linestrings are imported in the correct geometric coordinates but some other LINESTRINGs are imported with the negative longitude, so that LINESTRINGs are mirrored to the west.

When I open the DXF file with my QUANTUM GIS viewer I can see all the LINESTRINGs together with the correct coordinates, but when I export them to the postgres database and I print them in a map (by a web application) or later I export then I see some of the LINESTRINGs mirrored "away home" from the original coordinates.

That is my ogr2ogr commandline:

ogr2ogr -select "" -f PostgreSQL "PG:host=192.168.*.* user=**** password=**** dbname=****" C:\PATH\FILE.dxf -s_srs "EPSG:23031" -a_srs "EPSG:23031" -lco GEOMETRY_NAME=the_geom -nlt LINESTRING

There is no error output from the ogr2ogr tool, so it is difficult to get the reason of that problematic.

Any help or recommendation?

I have tried without much effect the ogr2ogr opions: -lco GEOM_TYPE=geometry -lco DIM=2 --config PG_USE_COPY YES -lco PRECISION=YES -clipsrc xmin ymin xmax ymax


texts:
I'm trying to generate a simple profile for some gridded topography data in GMT. However, even though the end of the segment produced by project is clearly in the NaN region of the data set (I usually mask it but included NaN as grey here). When I extract values with grdtrack it exports values for the entire length of the segment. The beginning of the segment also starts in NaN region successfully exports NaN until valid values are reached (second picture). I've  double checked that my projections and regions are consistent throughout.  Anyone know what might be causing this behaviour

gmt begin lake png
    gmt grdimage $DEM $region -JM6i -Qc
    gmt project -C-66.225/18.330 -A105 -L0/0.9 -G0.001 -Q | gmt grdtrack -G$DEM > $depth
    cut -f3,4 $depth > $depth1
    gmt plot $depth  -W0.1c,red
gmt end


texts: R 4.0.4 Raster package treating Albers as unknown projection

texts:
I'm running a fresh R installation with RStudio. Installed packages are:
library(randomForest)
library(randomForestExplainer)
library(raster)
library(rgdal)
library(sf)

When I try loading any raster with South America Albers Equal Area Conic projection (ESPG: 102033) I receive the following warning:

Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded datum Unknown based on Australian Natl & S. Amer. 1969 ellipsoid in Proj4 definition

Keeping the projection is important to me, cause I need to be able to collect point samples from that raster and I need to make sure the points and rasters are projected in the same coordinate system. Albers seems to be the most accurate for my study area.
I have packages updated to the last version availble @ cran and I cannot seem to find a solution for this.
Projection details are (quoting QGIS 3.4.4 defs):

Extent: Extent not known
Proj4: +proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs

I need a way to import rasters with albers projections without having them discarded.


texts: standard for generating a geolocation error ellipse from lines of position/bearing

texts:
Given SEVERAL points, where each point is a line of bearing (LOB)/line of position, and thus each point consists of a (1) latitude & (2) longitude position on the Earth's globe that this LOB is received from and (3) an angle of arrival in degrees from true north measured clockwise that this LOB was received from, as well as (4) a fixed length of each of these LOBs in nautical miles

From these SEVERAL points, I want to compute ONE location error ellipse consisting of one (1) latitude & (2) longitude coordinate which is the center of the error ellipse, (3) a semimajor and (4) semiminor distance in nautical miles of the ellipse, and (5) a orientation angle in degrees from true north to the semimajor axis of the ellipse.

I've been trying to look for the most commonly accepted/most popular method to do this, since this has definitely been done several times before.

I've seen this link: Standard Deviational Ellipse with Open Source Python (GDAL/OGR etc) and read a lot about the standard deviational ellipse, but I have a qualm with it.

If I were to use that standard deviation ellipse, I would probably compute all the points of intersection between each LOB, and then take all the intersection points (lat, long) and use those as (x, y) values as input to those equations to generate that ellipse.

My problem with this is that would assume that our (lat, long) pairs are in Cartesian space - which they aren't. If I were to assume so anyway (which I believe is invalid), then the resulting standard deviational ellipse would be in cartesian space, yet converting that to the earth's geocoordinates polar space wouldn't make sense to me.

If I were to NOT assume that they are in cartesian space, then I would convert the (lat, long) to (x, y, z), and the use this (x, y, z) to compute an error ellipsoid using those equations, then project the resultant ellipsoid on the earth's ellipsoid to get back to (lat, long) geocoordinate polar space, giving me a 2D (curved) error ellipse on the surface of the Earth. I don't have an altitude at each point necessary to compute the (x, y, z), but theoretically I can just use sea level so when I project the (x, y, z) ellipsoid back into (lat, long) geocoordinates, it can be altitude-agnostic. I can see this method making sense to me, but it seems a little too elongated for something that I believe is already solved.

I've looked at other papers as well, specifically those that generate an error ellipse just from two (a pair of) LOBs:
https://www.ngs.noaa.gov/PUBS_LIB/AlgorithmsForConfidenceCirclesAndEllipses_TR_NOS107_CGS3.pdf
and then papers that combine error ellipses together so I can get a final single error ellipse:
https://calhoun.nps.edu/bitstream/handle/10945/32273/96Sep_Orechovesky.pdf;sequence=1
but I've only scratched the surface on these in favor of looking for a less complex solution.

TL;DR what generally is the standard in this field to generate a geocoordinate error ellipse from several lines of bearing/position?


texts:
I'm working with ArcGIS 10.2. I used this code to select one record of an attribute table:

import arcpy
# Set the workspace
env.workspace = "F:ArcGIS\Trab\Model\Lin_150_tr.shp"
# Make a layer from the feature class
arcpy.MakeFeatureLayer_management("Lin_150_tr", "lyr5")

arcpy.SelectLayerByAttribute_management("lyr5", "NEW_SELECTION", ' "FID" = 2 ')

And it works!
Now I need to select no one record but some records of the same table, so I used a list of values and a loop:

import arcpy
env.workspace = "F:ArcGIS\Trab\Model\Lin_150_tr.shp"
arcpy.MakeFeatureLayer_management("Lin_150_tr", "lyr5")

crit_150=[2,10,200,323,349,384,433]
for a in crit_150:
    arcpy.SelectLayerByAttribute_management("lyr5", "NEW_SELECTION", ' "FID" = a ')

But it gives me an error:

Runtime error  Traceback (most recent call last):   File "<string>", line 21, in <module>   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 6494, in SelectLayerByAttribute     raise e ExecuteError: ERROR 000358: Invalid expression Failed to execute (SelectLayerByAttribute).

How can I do this?


texts: Change column to row using awk inside of GRASS command

texts:
I have 2000+ raster files and counting their statistics.
For now I am exporting all to Excel, but to work on it, I need better segregation.

My output:

Name1
1  123
2  243
3  455
Name2
1  432
3  655
4  234
Name3
2  344
3  345...

I have 4 possibilities: 1, 2, 3 and 4, but not all rasters have all combinations. I would like to use awk to create output like that:

name   1    2    3    4
Name1  123  243  455
Name2  432       655  234
Name3       344  345

So I was trying with this:

for red in {43..64}
do
for nir in {65..193}
do
r.stats -cnl 'gr_r'$red'_n'$nir | \
awk '$1 == 1 {print($2)}' ...

And how to create different rows?


texts: Converting u and v from WGS84 to NZTM2000 - Do I need to do vector conversion?

texts:
I have u and v wind component values in latitude and longitude in WGS84. I want to convert this to NZTM2000. I've converted the lat/lon in ArcView - (and will try in R) - but do I need to also do some sort of adjustment to the u and v values to account for the new coordinate system?

If so, is there a formula or perhaps a method in Matlab, ArcView, or R?


texts:
I have a sample of data taken from regions of a country! i want to test if there are any spatial autocorrelation on my data using Moran Indice test.
the null hypothesis: is that there is no spatial autocorrelation
here a sample of my data:

i used this steps to calculate it on r:

 library(RODBC)
 setwd('e:/r/moran')
 channel <- odbcConnectExcel('moran.xls')
 data <- sqlFetch(channel, 'wilaya')
 inf.dists <- as.matrix(dist(cbind(lon=data$Lon, lat=data$Lat)))
 inf.dists.inv <- 1/inf.dists
 diag(inf.dists.inv) <- 0
 library(ape)
 Moran.I(data$year2009, inf.dists.inv)

and i got these results:

$observed
 -0.02229578

$expected
-0.02702703

$sd
 0.03455708

$p.value
 0.8911011

the $observed ~= $expected and they are negative! does that mean there is a little dispersion! but we cannot reject the null hypothesis, so there is no spatial autocoorelation between the regions!

the problem is that i'm not sure about my interpretation ?
and i want to know if the method i used for calculating the weight matrix for moran.I function is right?


texts: Compiling Python scripts (to .exe) that use ArcGIS Geoprocessing Tools?

texts:
I've been coding with Python for several months now and have developed some reasonably complex scripts for primarily geoprocessing tasks. That being said, I'm still learning a lot as I'm coming from a SQL/VBA/VBScript background.

I know that compiled code typically runs faster than code that must be processed by a language interpreter, so I'm interested in the possibility of compiling a geoprocessing Python script to a .EXE file for working with big data.

Is this even possible? If it is, what is the best way to compile a Python (.py) script that is importing the arcgisscripting or arcpy modules?

I spent a few minutes trying to find what I want to do and the search returned this article among others: http://www.ehow.com/how_2091641_compile-python-code.html

The compiler seemed to work, but upon executing the resulting .EXE file, it gave a cryptic error saying some files were unavailable.

The Python script runs what seems to be reasonably well from the command line, but I'm wondering if I could see some slight improvement if I were able to compile the .py file. Again, I'm working with some big datasets that are taking +20 hours to process (delineating watersheds from input water-quality sample sites). I'll take anything I can get in way of improvements.

The script ran 10% quicker outside of ArcGIS from the command line using a test set of sites versus setting the script up as a script tool in a new toolbox in ArcCatalog. I've been running the script from the command line w/o any instance of ArcGIS open on a dedicated machine.

So, is it possible to compile Python scripts that import the arcgisscripting module and that call ArcToolBox tools?

EDIT

Thanks for the input, this is helpful for me. The script is largely a way to coordinate a number of ArcGIS tools and to output in desired formats/locations/with appropriate attribution. I've already trimmed some fat I think by writing to a scratch folder instead of a scratch personal geodatabase for some interim raster files so they can be stored in the ESRI GRID format vs. the IMG format. I'll check out the profiler suggestions though.

There are some in my office that question Python saying "that compiled code is so much quicker than code running through an interpreter" mainly in comparison to, say, a compiled Visual Basic program or VB.NET program, but that is a good point that the tools are going to take time either way. And, it seems like with present day computing machines that interpreted code may not be that much slower than compiled code to warrant going that extra mile.

EDIT - update on optimization of the program with raster formats.

Wanted to follow up on my "optimization" of this Python program, and I was able to shave 2 hours of processing time by writing interim rasters to GRID format instead of to a personal geodatabase. Not only that, there was a SIGNIFICANT reduction in data size disk space consumption. The original run I did writing all rasters (and they were only point features converted to rasters, and then watershed rasters) resulted in 37.1 GB of data just for those files. Writing the latter two data outputs to a folder in GRID format was reduced to 667 MB of data.

I'd be curious to see how a file GDB would handle these data though mainly in way of the size of the data. But, cutting my processing time down from 9.5 hours to 7.5 hours certainly is enough to advocate for dealing with rasters outside of geodatabases in the GRID format.


texts:
This question is relative to a past question of mine at Rendering raster layer from DTED based on height color map (that remained unanswered), so if you believe that I should merge those two, let me know.

I have a MADtedLayer created from a dted catalog using military analyst. The only way I managed to display the data, was using the UseUniformRenderer property of the layer. I want to be able to change the renderer.
To be more precise, I want to display each elevation with a specific color, that is already available from a file.

My guess is that I need a RasterDiscreteColorRenderer with a corresponding RasterColormap (although I'been unsuccesful with a StretchRenderer also)

Here I provide my un-succesful so far code:

public IRasterRenderer CreateDtedRenderer()
{
    IRasterColormap colormap = MakeDefaultColormap();

    IRasterDiscreteColorRenderer renderer = new RasterDiscreteColorRenderer();
    renderer.Colormap = colormap;
    return (IRasterRenderer) renderer;
}

public IRasterColormap MakeDefaultColormap()
{
    var clrString = DtmColorTable.GetClrAsString();
    var blues = new List<double>();
    var reds = new List<double>();
    var greens = new List<double>();
    using (var reader = new StringReader(clrString))
    {
        string line;
        while ((line = reader.ReadLine()) != null)
        {
            var vals = line.Split(' ');
            reds.Add(Int32.Parse(vals[1]) / 255.0);
            blues.Add(Int32.Parse(vals[2]) / 255.0);
            greens.Add(Int32.Parse(vals[3]) / 255.0);
        }
    }

    IRasterColormap colormap = new RasterColormap();
    colormap.BlueValues = blues.ToArray();
    colormap.RedValues = reds.ToArray();
    colormap.GreenValues = greens.ToArray();
    return colormap;
}

public static IMADtedLayer CreateMADtedLayer(IWorkspace workspace, string catalogName)
{
     var inputDatasetName = System.IO.Path.GetFileName(catalogName);

     // Open the raster catalog
     var rasterCatalog = ((IRasterWorkspaceEx) workspace).OpenRasterCatalog(inputDatasetName);

     // Create the MADtedLayer
     IMADtedLayer maDtedLayer = new MADtedLayer();
     if (!maDtedLayer.Setup((ITable) rasterCatalog)) return null;

     ((IGdbRasterCatalogLayer) maDtedLayer).Renderers.Insert(0, CreateDtedRenderer());
     maDtedLayer.UseUniformRenderer = true;
     maDtedLayer.DisplayElevation = true;
     maDtedLayer.ResamplingType = rstResamplingTypes.RSP_NearestNeighbor;
     ((IWireframeDisplayProps) maDtedLayer).DisplayWireframe =
            esriMAWireframeDisplayEnum.esriWireframeDisplayNone;

     return maDtedLayer;
}

Has anyone ever changed the rendere of a MADtedLayer (or a IGdbRasterCatalogLayer in general?)


texts: Cache on demand feature for cached map services

texts:
I have a cached map service that I am using in a Silverlight based application. For this cached service, I have cached a core area and checked the option for On Demand Cache. I find that the areas within On Demand Cache are very slow. I am trying to understand how this process works. If I am zoomed in at the scale of 4000, it should be generating cache for a very small area, not sure why it takes much longer. I am using ArcGIS online tiling scheme. Does it generate the cache for all scale levels or just the one level.

I am looking for some feedback.

Thanks
Jay


texts: Simultaneous raster downsampling/downscaling and upsampling/upscaling for false accuracy minimisation

texts:
I want to have a raster stack of two different resolutions.
One is from VSNIR and one from SWIR sensors (hyperspectral).
The VSNIR has a resolution of 25 cm and the SWIR of 80 cm and both data have the same importance for my study.
If I upscale the VSNIR and downscale the SWIR rasters to 60 cm, the resulting stack could have the minimal loss of accuracy from both bands. The number 60 comes as the area of the square which is the average of the areas of the two pixels 0.5 x √(25 x 25)+(80 x 80). I have doubts, that I have to consider also that the upscaled raster will have a greater loss of accuracy (the result is extrapolated) and add a parameter which changes the optimal pixel resolution.

Any ideas of what would be the optimal common resolution for bands of such different resolutions? I will use this data for spectral regression analysis.


texts:
I opened a sub-dataset of an HDF file (NDVI layer of MODIS HDF file) in QGIS, but pixel values are beyond the valid range (-2000,10000). I see many pixels way less than -2000 and way over than 10000. There is something wrong here.

When I read layer into numpy array in Python, all pixel values between valid range, so everything is expected.

you can download the file here, MODIS hdf file Extract the first layer, NDVI, no=0.

information about QGIS version on Ubuntu 12.04 LTS:
QGIS version    2.8.1-Wien  QGIS code revision  exported
Compiled against Qt 4.8.1   Running against Qt  4.8.1
Compiled against GDAL/OGR   1.10.0  Running against GDAL/OGR    1.10.0
Compiled against GEOS   3.4.2-CAPI-1.8.2    Running against GEOS    3.4.2-CAPI-1.8.2 r3921
PostgreSQL Client Version   9.1.15  SpatiaLite Version  4.1.1
QWT Version 5.2.2   PROJ.4 Version  480
QScintilla2 Version 2.6.1


texts: GDAL translate -AttributeError: 'NoneType' object has no attribute 'RasterCount'

texts:
I am trying to convert a netCDF file to geotiff, I am getting he following error at line 14. I just can't seem to figure out what the problem is, I will put the whole code below. Please help,

AttributeError                            Traceback (most recent call last)
<ipython-input-14-d0e422b9cf8e> in <module>()
     12     driver = gdal.GetDriverByName(fileFormat)
     13     band = 1
---> 14     numBands = dataset.RasterCount
     15
     16     # Write gdal_translate into multiple bands - again specific to my .nc file

AttributeError: 'NoneType' object has no attribute 'RasterCount'

The code :

Step 1: GDAL Translate to convert to TIFF

import gdal
import datetime

for i in range(1980,2010):
    #change this to suit your .nc filename pattern
    fileName = 'full_data_v6_precip_05' +str(i)+ '.nc'
    year = i
    startDate = datetime.date(year, 1, 1)
    # create new file to write command prompt lines to
    file = open('commandpromptoutput.txt','a')

    #parse .nc file to count how many bands
    fileFormat = "NETCDF"
    dataset = gdal.Open(fileName,gdal.GA_ReadOnly)
    driver = gdal.GetDriverByName(fileFormat)
    band = 1
    numBands = dataset.RasterCount

# Write gdal_translate into multiple bands - again specific to my .nc file
    for i in range(numBands):
        file.write("gdal_translate -a_srs EPSG:4326 -b "+str(band)+" "+ fileName +" "+ str(startDate) + ".tif\n")
        startDate = startDate + datetime.timedelta(days=1)
        band += 1

    file.close()

Step 3: Upload to EE from GCloud

from datetime import datetime
import os
path = "C:\WorkingDirectory"
# create new file to write command prompt lines to
file = open('commandpromptoutput.txt','a')
epoch = datetime.utcfromtimestamp(0)
def unix_time_millis(dt):
    return (dt - epoch).total_seconds() * 1000.0

for root, dir, files in os.walk(path):
    for filename in files:
        filedate = datetime.strptime(filename[0:-4], '%Y-%m-%d')
        filetime = unix_time_millis(filedate)
        file.write('earthengine upload image --asset_id users/xxx'+str(filename[0:-4])+' --nodata_value -999 --property system:time_start='+ str(int(filetime))+ ' gs://xxx' + str(filename)+"\n")

file.close()


texts: Recovering toolbar from disconnected window?

texts: I've faced this problem many times. I created a python addin toolbar and installed it in ArcGIS 10.2 . I work on dual monitors in the office. If accidentally this toolbox is left on the secondary screen and I open my laptop again and disconnect the second screen, I can't see the toolbox, it is hidden or remains on second monitor. I can't find a way to move this toolbox into primary screen. Can any body help me with this weird situation ?


texts: Convert Google Maps KML to preferred format on a Hyundai Ioniq EV's BlueLink Map

texts:
I want to upload data for all Eduroam WiFi locations in Canada to my car navigation system which has an "Import Address Book from USB device" option.
I downloaded the current locations I have on my car and it looks like this:
<?xml version="1.0" encoding="utf-8" ?>
<kml>
    <Folder>
        <Placemark>
            <name>Eduroam</name>
            <Point>
                <coordinates>-123.143989,49.179125</coordinates>
            </Point>
        </Placemark>
        <Placemark>
            <name>Edroam2</name>
            <Point>
                <coordinates>-123.152900,49.155589</coordinates>
            </Point>
        </Placemark>
        <Placemark>
            <name>Home</name>
            <Point>
                <coordinates>-123.130431,49.191978</coordinates>
            </Point>
        </Placemark>
    </Folder>
</kml>

So I think it wants latitude and longitude.
Where as the custom Google map provided by Eduroam Canada looks like:
<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <Document>
    <name>EduroamCanada</name>
    <description/>
    <Style id="icon-22-normal">
      <IconStyle>
        <scale>1.1</scale>
        <Icon>
          <href>https://www.gstatic.com/mapspro/images/stock/22-blue-dot.png</href>
        </Icon>
        <hotSpot x="16" xunits="pixels" y="32" yunits="insetPixels"/>
      </IconStyle>
      <LabelStyle>
        <scale>0</scale>
      </LabelStyle>
    </Style>
    <Style id="icon-22-highlight">
      <IconStyle>
        <scale>1.1</scale>
        <Icon>
          <href>https://www.gstatic.com/mapspro/images/stock/22-blue-dot.png</href>
        </Icon>
        <hotSpot x="16" xunits="pixels" y="32" yunits="insetPixels"/>
      </IconStyle>
      <LabelStyle>
        <scale>1.1</scale>
      </LabelStyle>
    </Style>
    <StyleMap id="icon-22">
      <Pair>
        <key>normal</key>
        <styleUrl>#icon-22-normal</styleUrl>
      </Pair>
      <Pair>
        <key>highlight</key>
        <styleUrl>#icon-22-highlight</styleUrl>
      </Pair>
    </StyleMap>
    <Style id="icon-22-nodesc-normal">
      <IconStyle>
        <scale>1.1</scale>
        <Icon>
          <href>https://www.gstatic.com/mapspro/images/stock/22-blue-dot.png</href>
        </Icon>
        <hotSpot x="16" xunits="pixels" y="32" yunits="insetPixels"/>
      </IconStyle>
      <LabelStyle>
        <scale>0</scale>
      </LabelStyle>
...

You can find the eduroam Canada custom  map here: https://www.google.com/maps/d/u/0/viewer?mid=1mgS5P6gcv0A_5-QuKUb3pNg9VX8&ll=49.1560847597644%2C-123.14908251417819&z=13
How can I convert the Google Map format to what Bluelink wants?
Right now it just says "No supported files" when I stick the USB in.


texts: Extracting a region's outline of roads from a GIS database, as a vector image

texts:
It's now about 5 days I've been trying to extract an outline of the internal transportation routes of a city (Tehran) from various kinds of GIS sources, and exporting that outline as a vector image, preferrably as an SVG.
To that end, I have tried using QGIS for opening the tileset in https://openmaptiles.com/downloads/dataset/osm/asia/iran/tehran/#8.91/35.102/51.0771 and every single GIS file at https://download.bbbike.org/osm/bbbike/Tehran/ . I have ransacked the Internet for help, but so far I've been totally clueless as to where to even begin.
I have attached the result I'm trying to achieve at the end of this post.


texts: How to compile QGIS with Oracle Spatial support for Linux?

texts:
I compiled it on Fedora and Ubuntu, but got error "Cannot get Oracle select dialog from provider." while trying to add Oracle layer.
Packages in repositories don`t have oracle spatial support.
Can somebody help me?
I need deb or rpm package with oracle spatial and qgis-mapserver support.


texts: Using DEM to find permanently shadowed craters

texts:
I'm using ArcGIS Pro to try to identify permanently shadowed craters on the North Pole of the Moon. Pictured below is a screenshot of the DEM I'm using which was collected from the laser altimeter on the Lunar Reconnaissance Orbiter. It is available for download from USGS here:

I can use the Hillshade tool with a low altitude angle and model_shadows turned on to create a raster like the one below that shows the shadowed regions of the craters for those particular illumination conditions. The picture below shows the results of a Hillshade with Azimuth 315 deg and Altitude 5 deg.

In a simplified world, I could just run this hillshade for a variety of Azimuth inputs and sum the resulting rasters. The pixels that still had a value of 0 would be permanently shadowed regions. There is one problem with this plan though. The moon is not flat, especially over such a large region. This DEM covers the whole surface of the Moon north of 75 deg latitude and even lower at the corners. Because of this, the elevation of the sun would vary quite a lot for various pixels in the image.

To summarize:
Is there a way to calculate shadows caused by terrain on a large, spherical surface?

The Hillshade tool operates only on planar maps.

EDIT: Esri has a help page where they basically give you the hillshade formula and explain all the terms. I was able to modify this formula by calculating the unique azimuth and zenith angle for each pixel. The result below is almost what I'm looking for. In this image, the sun is off the top of the page and it is easy to see how the curvature of the moon makes the other side darker due to the higher zenith angles. The only thing this image is missing is shadows, which is why some of the crater walls on the bottom image are still illuminated despite being on the far side from the light. Unfortunately, I can't find any documentation on how shadows are calculated.