texts: Selecting the correct projection for NAD_1983_HARN_WISCRS_Bayfield_County_Feet

texts:
I am trying to map the following projection but having issues finding a one that works. I ahve used http://prj2epsg.org/search which gives a few results but none seem to work that well.

PROJCS["NAD_1983_HARN_WISCRS_Bayfield_County_Feet",
GEOGCS["GCS_North_American_1983_HARN",
DATUM["D_North_American_1983_HARN",
SPHEROID["GRS_1980", 6378137.0, 298.257222101]],
PRIMEM["Greenwich", 0.0],
UNIT["degree", 0.017453292519943295],
AXIS["Longitude", EAST],
AXIS["Latitude", NORTH]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["central_meridian", -91.15277777777777],
PARAMETER["latitude_of_origin", 46.66964837722222],
PARAMETER["scale_factor", 1.0000331195],
PARAMETER["false_easting", 750000.0010000002],
PARAMETER["false_northing", 487372.65900000004],
UNIT["foot_survey_us", 0.3048006096012192],
AXIS["X", EAST],
AXIS["Y", NORTH]]

Use Prj2EPSG with the text above (which is what geoserver sees as the native SRS), I get matched with 2236,2241,2239,2254,2257,2260..etc

When I use the prj file instead, I get this.

PROJCS["NAD_1983_HARN_WISCRS_Bayfield_County_Feet",
GEOGCS["GCS_North_American_1983_HARN",
DATUM["D_North_American_1983_HARN",
SPHEROID["GRS_1980",6378137.0,298.257222101]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic"],
PARAMETER["False_Easting",750000.001],
PARAMETER["False_Northing",487372.659],
PARAMETER["Central_Meridian",-91.15277777777779],
PARAMETER["Standard_Parallel_1",46.66964837722222],
PARAMETER["Scale_Factor",1.0000331195],
PARAMETER["Latitude_Of_Origin",46.66964837722222],
UNIT["Foot_US",0.3048006096012192]]

These match with 3857,24370..etc Obviously the two are different. Geoserver must be parsing it a bit different than prj2epsg does. prj2epsg uses the exact text of the prj file.

After trying to set a few of the declared SRS. When I go to Data->Layer Preview. The layer will show up how it should look just fine. However, when I go to Tile caching -> Tile Layers, both the 4326 and 900913 don't show anything. Typically they show up just fine with most of my other layers. Granted those were easier to use since they were 4326 to start with.

Any suggestions would be helpful.

EDIT: When I select 2236 for example. I can get it to work. However, it shows up in flordia instead of wisconsin. That makes sense since it is pretty much trying to project it on the map in florida. Sounds like I need a custom projection for this? Not sure?


texts: Search a certain value in a cursor in multiple fields

texts:
I am trying to find the phrase REQUIRE!!! in certain fields (PROJ_NAME, MAJOR_LITH, OWNER_ID) in a feature class. If it has the value REQUIRE!!! it will send user an email. If there are no value that is REQUIRE!!! it will update another feature class.

I got the updating part but could not figure out the searching field part. I can't figure out what is wrong with the cursor part.

    with arcpy.da.SearchCursor(Strip_Geology_BOB6, ['PROJ_NAME', 'MAJOR_LITH', 'OWNER_ID']) as cursor:
for row in cursor:
    value = ['REQUIRED!!!']
    if row in value:
        MESSAGE_REGION = "From: %s\r\n" % FROM + "To: %s\r\n" % TO + "Subject: %s\r\n" % "Strip Geology - Missing Required Fields" + "\r\n" + "Auto generated Message.\n\rPlease fix required fields in Strip Geology Feature Class"
        server = smtplib.SMTP(SERVER)
        server.sendmail(FROM, TO, MESSAGE_REGION)
        server.quit()
    else:
        arcpy.Select_analysis (Strip_Geology_BOB6, Strip_Geology_BOB6_Select, "DATAUPLOAD = 'No'")
        arcpy.Append_management(Strip_Geology_BOB6_Select, SDE_Strip_Geology, "TEST", "", "")
        arcpy.CalculateField_management(Strip_Geology_BOB6, "DATAUPLOAD", "\"Yes\"", "PYTHON", "")
        arcpy.Delete_management (Strip_Geology_BOB6_Select)


texts:
im working with with QGIS 2.14.2 server version with postgres 9.5 database.
i want to know how to import SVG symbols stored in QGIS library  into POSTGRES?


texts: Classification based on NDVI values of two images

texts:
I want to find the areas where the NDVI value were higher then 0.4 at the end of April and NDVI values were lower then 0.3 at the end of June. I have tried to make it using below demonstrated script, unfortunately it did not work out. Could you please assist me on improving the script below or suggest me another ways to do it.

var test = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR/")
                .filterBounds(point)
                  .filterDate('2000-04-15', '2000-04-30')
                    .normalizedDifference(['B4', 'B3'])
                      .gt(0.4). and (
            ee.ImageCollection("LANDSAT/LT05/C01/T1_SR/")
                .filterBound(point)
                  .filterDate('2000-07-15', '2000-07-30')
                    .normalizedDifference(['B4', 'B3'])
                      .lt(0.3));
print(test);
Map.addLayer(test);


texts: Creating heat map weighted by population in ArcGIS Online?

texts:
I want to create a hotspot map of incidents by City but I would like to weight it by population.

Is this possible in ArcGIS Online?


texts: Filter table of attributes by numerical values (all values between x and y)

texts:
I am trying to select features based on their area (in the column "surface") in the table of attributes.
I wrote: SELECT "surface" BETWEEN 300 AND 1500, it doesn't work.
I also tried: "surface" >300 <1500 but it doesn't work either.


texts: Export geography field from SQL for use in GIS

texts:
I currently have a table in SQL server management studio (2012) with a geography column.

It definitely has valid spatial data, as I can view it in SSMS:

What I want is to export the data as a flatfile csv, so that I can load it into GIS software, specifically QGIS. I tried exporting as csv and loading the csv into QGIS, selecting the geography column as WKT, however QGIS did not recognise it. Is it possible to convert whatever format SQL stores Geography information into WKT?


texts: GeoServer layer preview doesn't work for chinese columns or columns that has underscore in name

texts:
I am learning to use GeoServer. I published a layer from a shp, which I created by importing a csv file into QGIS in utf-8 encoding, and exporting as shp in utf-8 encoding. I added a style to the layer. I got the style from geoserver cookbook for point with default label. The only thing I changed is the contents in <ogc:PropertyName>, as I have to set it to the name of the column I want to show. I changed it into school_name. But then the layer doesn't show in layer preview.

Things I've tried:

change <ogc:PropertyName> content to the column lat. This column contains latitude. The type shown on GeoServer is double. This works.
change <ogc:PropertyName> content to school. This column contains Chinese characters. The type shown on GeoServer is string. The layer shows but the labels are ??.
change <ogc:PropertyName> content to school and change the xml encoding to utf-8 and UTF-8. Again, the layer shows, but the labels are ??
change <ogc:PropertyName> content to school and add format_options=CHARSET:UTF-8 in the WMS request. Again, the layer shows, but the labels are ??.
change content to postal_code. This column contains postcodes which are just numbers. The type shown on GeoServer is long. The layer doesn't show.

So how do I show Chinese labels? And how do I fix the problem where column names with underscore doesn't show?


texts: Is there anyway to read grib file without download it?(directly from the link)

texts:
grbs = pygrib.open(
   'https://www.ncei.noaa.gov/data/north-american-mesoscale-model/access/historical/forecast/202001/20200101/nam_218_20200101_0600_024.grb2')
grb = grbs.read(1)[0]

This is what I am trying to do, but it gives me the following error
The link has no issue, since i can use that link and download the file in the local machine, and use pygrib.open to read file from local machine, But now I am having trouble to read file directly from the link. I read through pygrib documentation, it seems like they did not provide a way to do that. Is there any way can read grib file drectly from link?


texts: Get WMS service properties for cql filtering?

texts:
From a wms service I know its url and a layer name called: CATAST_Pol_Municipio

https://idena.navarra.es/ogc/wms

I need to apply a cql_filter on a GetMap request on that layer, filtering by national cadastral reference. Tried cql_filter=nationalCadastralReference=xxxxxxxx but it says that there is no property with nationalCadastralReference name.

  Rendering process failed PropertyDescriptor is null - did you request
  a property that does not exist?

How can I know the wms properties that I can filter the request with?


texts:
When I view MOD02 data in ENVI directly with ...
File--Open image file, ENVI automatically loads data.
So what is the meaning of Emissive?

The description of MOD02 data can be found in MOD 02 - Level-1B Calibrated Geolocation Data Set and MODIS Data Product Non-Technical Description - MOD 02, 03.


texts: Layers from GRASS db won't load in QGIS 2.18.4 desktop

texts:
I have been learning QGIS online but now I am stuck at a module. I am supposed to launch a project that includes GRASS vector and point layers. However, when I do so I always receive the 'Handle bad layers box'. When I use the browser panel and navigate to the the layers, they appear with red 'X' icons.

The same thing happens with the QGIS Alaska GRASS demo. I have GRASS 7 plugin enabled.

Is there something I am missing?


texts: Projected coordinate reference system for aviation navigation and instrument approach charts

texts:
I am trying to find out if aviation organizations like the ICAO have written guides for what projection CRS must be used in the following charts:

High/low altitude navigation charts
IFR approach procedures charts

I am using QGIS.


texts: Copy features from a SQLite multipolygon to SQLite polygon layer

texts: I need to copy features from a multipolygon layer that contains multipart features into a polygon layer in SQLite in QGIS. I know I can run the "Multiparts to Singleparts" algorithm and simply copy and paste but this is something I'm going to have to do a lot and was wondering if there is some easier way.


texts:
I've found several examples for how to draw a simple arrow (on one line end). Also, I've found few that explain how it should work to put one on the other end of a line. Unfortunately I don't have any success with that.

The code looks as follows:

        <Rule>
          <Title>line with arrow on both endings</Title>
          <LineSymbolizer>
            <Stroke>
              <CssParameter name="stroke">#E60000</CssParameter>
              <CssParameter name="stroke-width">1</CssParameter>
            </Stroke>
          </LineSymbolizer>
          <PointSymbolizer>
            <Geometry>
              <ogc:Function name="startPoint">
                <ogc:PropertyName>the_geom</ogc:PropertyName>
              </ogc:Function>
            </Geometry>
            <Graphic>
              <Mark>
                <WellKnownName>shape://oarrow</WellKnownName>
                <Fill>
                  <CssParameter name="fill">#E60000</CssParameter>
                  <CssParameter name="fill-opacity">1.0</CssParameter>
                </Fill>
              </Mark>
              <Size>30</Size>
              <Rotation>
                <ogc:Add>
                  <ogc:Function name="startAngle">
                    <ogc:PropertyName>the_geom</ogc:PropertyName>
                  </ogc:Function>
                  <ogc:Literal>180.0</ogc:Literal>
                </ogc:Add>
              </Rotation>
            </Graphic>
          </PointSymbolizer>
          <PointSymbolizer>
            <Geometry>
              <ogc:Function name="endPoint">
                <ogc:PropertyName>the_geom</ogc:PropertyName>
              </ogc:Function>
            </Geometry>
            <Graphic>
              <Mark>
                <WellKnownName>shape://oarrow</WellKnownName>
                <Fill>
                  <CssParameter name="fill">#00E600</CssParameter>
                  <CssParameter name="fill-opacity">1.0</CssParameter>
                </Fill>
              </Mark>
              <Size>30</Size>
              <Rotation>
                <ogc:Add>
                  <ogc:Function name="endAngle">
                    <ogc:PropertyName>GWAI_GEOM</ogc:PropertyName>
                  </ogc:Function>
                  <ogc:Literal>0.0</ogc:Literal>
                </ogc:Add>
              </Rotation>
            </Graphic>
          </PointSymbolizer>
        </Rule>

What I get is this:

Two lines with start- and end-arrow both on the same line end and end-arrow not turned at all. What is the error I have in my code?

GeoServer Version: 2.13.2

/EDIT: When removing the "Rotation" bit from endPoint Graphic, the result looks the same.


texts:
Using ArcMap 10.1, I’m editing values in a feature class which was sent to me from within my organization. I realized I had a problem when trying to change the symbology of this feature class- the values I updated appear in the attribute table but some don’t appear in the Properties dialog box when I try to change the Symbology to categories/unique values (trying to visually highlight the new values I've added).

The structure of the attribute table cannot be changed as it was created to link with other systems. The problem is only occurring with one field; other fields which I updated do have their unique values appear in the Properties dialog box. I checked to see if this was created to be a Read-only field and it was not... Is there some protection on the feature class I don't know about? Did the edits I made somehow not get stored? Edits were made using field calculator + some manual edits in a Edit session which have all been saved.


texts: Trying to update field in a feature Without "attInspector" ...Using ArcGis JavScript api

texts:
I cannot find any changes in the field. Please help me out.

I am attaching my code:

for (var i = 0; i < features.length; i++)
{
  feature = features[i];
  feature.attributes["X"] = Value;

  updateFeature.getLayer().applyEdits(null, [feature], null);
}


texts: Can arcpy read dynamic text as appeared in ArcMap?

texts:
Just wondering if we can read dynamic text in arcpy?
For example, assume we have a Text Element called myText which contains a piece of dynamic text Map <dyn type="page" property="index"/>. Reading this Text Element is as below:
text = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "myText")[0].text

text is now a string of Map <dyn type="page" property="index"/> which is not what I actually want. In ArcMap, this Text Element is showed as Map 1, and this is the value I want for my text variable.
Is there any simple way to achieve this with arcpy?
As suggested in Access resolved/rendered dynamic text from ArcPy or ArcObjects?, this can be done by either parsing the text (which is not easy to do properly) or saving the mxd to v9.3 (which would slow down the process and may pose a compatibility risk)
UPDATE 1
Decided to go ahead with parsing the text as the dynamic texts in my MXDs are not that complicated. Cannot afford the approach of saving the mxd to v9.3 due to the huge number of MXDs to deal with.
Anyway, any other suggestions are more than welcome.


texts: Changing draw text element to be right to left using ArcPy?

texts:
I'm trying to produce a set of maps in certain locations using Data Driven Pages. I'm using a script written in python that generates a dynamic table in each sheet. The table has attributes of all the entities participating in that particular sheet.

The dynamic table is been generated just fine expect I want the attribute to be written in right to left writing instead of default left to right.

How to generate text element that is written in right to left?

The function needs to be as in the label's editor option:.

This question was asked in different variation as Understanding UnicodeDecodeError: 'utf8' codec in ArcPy script? here and pretty much the same in GeoNet. It didn't get any useful answers.


texts: Why is ogrinfo showing no layers from an ESRI FileGDB?

texts:
I have an ESRI File Geodatabase (FileGDB) (size is 15+GB), which has over 10 layers, when I open it up in ArcGIS, everything seems to be working.

I want to move this database to postgreSQL with OGR ogr2ogr.

First I am checking the ogrinfo with:

ogrinfo pathToFileGDB

I get result:

INFO: Open of `C:\Foo\My Documents\ArcGIS\Default.gdb'
  using driver `FileGDB' successful.

But I see no layers, nor does the ogr2ogr see any layers. I should see all the layers here like I saw in ArcGIS. So I am unable to move it to postgreSQL.

I tested the same thing on a smaller FileGDB(800MB) and I can see all the layers with ogrinfo, also the ogr2ogr works.

What could be the problem? Is my DB too big, if so, what could be a workaround?
Any other suggestions?


texts: How to Create Voronoi Polygon Diagram within surrounding polygon?

texts:
I want to create a Voronoi Diagram within an outer polygon.

I want to retrieve a result very similar to this solution, which uses ST_DelaunayTriangles() to create Triangles.

Actually, it would suffice to decompose the polygon into multiple polygon with the blue markers as intersecting points:

I can't think of any approach to get the desired polygon set.
The most crucial part is that the outer polygon has way much more points that describe the polygon than I need to create the polygon subsets.
Simply using ST_DelaunayTriangles() would lead to several more polygons in the subset than are required and would contain even triangles outside the surrounding polygon, what I want to avoid at all.

This image shows an example of the problem:

This is an illustration of the result I want to retrieve:

Apart that, the outer lines are more complex, as you can see in the first image of the yellow outer polygon.

Any ideas?

Software Environment:

Java Backend with JDBC Connection to PostGIS Database.
Frontend is a usual Web Application running with node.js.


texts: Enable Mosaic Native Acceleration in GeoServer 2.17

texts:
I am trying to optimise my server's GeoServer configuration for serving WMS. My server is running Ubuntu 19.10. I'm following the 'GeoServer on steroids' guide which recommends enabling 'Mosaic Native Acceleration', amongst other things.

As I understand it, to do this I need to switch to native JAI, rather than the pure java version. I downloaded and installed the relevant files (following this guide) and also removed the pure java implementations from GeoServer's WEB-INF/libs directory as they were causing a conflict that prevented GeoServer from starting (also mentioned in the above guide).

My Geoserver status is as follows:

Data directory  /usr/share/geoserver/data_dir
Locks   0
Free locks
Connections     0
Memory Usage    194 MB / 672 MB
Free memory
JVM Version     Private Build: 1.8.0_252 (OpenJDK 64-Bit Server VM)
Java Rendering Engine   sun.java2d.pisces.PiscesRenderingEngine
Available Fonts     GeoServer can access 314 different fonts. Full list of available fonts
Native JAI  true
Native JAI ImageIO  true
JAI Maximum Memory  336 MB
JAI Memory Usage    0 KB
Free memory
JAI Memory Threshold    75%
Number of JAI Tile Threads  7
JAI Tile Thread Priority    5
ThreadPoolExecutor Core Pool Size   5
ThreadPoolExecutor Max Pool Size    10
ThreadPoolExecutor Keep Alive Time (ms)     30000
Update Sequence     578

I'd hoped this would allow me to check the 'Mosaic Native Acceleration' tick box on the 'Image Processing' setting page but — alas — the check box remains greyed out. How can I get this working?


texts:
I am using shortest path point to point and I have found the following problem:
I am using a road network which includes some roads that have one way and others are bidirectional.  My origin and destination points are not in the network but they are near the network.
If I run the shortest path algorithm without specifying a road direction I get the shortest path between the points, but if I define the one way roads, I do not get a result.
The problem seems to be that the shortest path algorithm first looks for the shortest distance between the points and the road network and then runs the algorithm using that information.  If the shortest distance between the points and the road network lies on a one way road that is not reachable then it gives a not path found error, even if there is a reachable road within the network only a few meters away.
Is there a way to modify the algorithm to search for other points? (Something like, if no route is find, search for the second closest point between the road network and the origin destination points and try again)


texts:
I have some csv lookup tables which are too heterogeneous for ordinary join so I want to process them in PyQGIS. However the csv files have no geometry column whatsoever. If used as QgsVectorLayer the layer ends up empty:

uri='path/to/csv/file.csv?delimiter=;'
csvLyr = QgsVectorLayer(uri, "csvLayer", "delimitedtext")
QgsProject.instance().addMapLayer(csvLyr) #throws 'source data cannot be found'
csvLyr.getFeature(0) #returns a non callable QgsFeature

Is there a way in PyQGIS to query and manipulate non-spatial csv data just like any other vector data (respectively just like in the GUI) or am I required to use the csv module/another Python library?

I am using Python/PyQGIS 3 in QGIS 3.6. The csv files all have header rows.


texts:
Seems like such a simple plugin, but it just won't work.

I'm trying to run Zonal Statistics to find the mean value of a raster's pixels within the area of a vector polygon, but all it will give me is 'NULL'. The polygon and the raster are both in the same CRS. There is one null data point within the area, it has a no data value of -9 and is identified as "null (no data). But I don't believe that that is the issue because choosing an area that does not contain a null still doesn't work.

In order to see the new columns in the attribute table I find it necessary to toggle editing on the polygon on and off. I've tried running it with editing both on and off, one time it worked with editing on! But never again. Is it supposed to require editing to be on?

I'm using QGIS 1.8.0 on a Mac running OS 10.7.5.


texts: Simplifying and plotting polygons in Leaflet package in R?

texts:
I was finding the plotting of shapefiles very slow in R. After reading this (How to speed up the plotting of polygons in R?) I found all the tips were helpful for plotting in base R. The second one (creating a reduced shapefile by removing small polygons using a custom function) was particularly useful.

library(rgdal)
library(leaflet)
library(rgeos)

full.shapefile <- readOGR(dsn="170411_maps", layer="SA2_2011_AUST")

simplified.shapefile <- gSimplify(full.shapefile, tol=0.01, topologyPreserve=TRUE)
simplified.shapefile = SpatialPolygonsDataFrame(simplified.shapefile, data=full.shapefile@data)

# Remove the islands
getSmallPolys <- function(poly, minarea=0.01) {
  # Get the areas
  areas <- lapply(poly@polygons,
              function(x) sapply(x@Polygons, function(y) y@area))

  # Quick summary of the areas
  print(quantile(unlist(areas)))

  # Which are the big polygons?
  bigpolys <- lapply(areas, function(x) which(x > minarea))
  length(unlist(bigpolys))

  # Get only the big polygons
  for(i in 1:length(bigpolys)){
    if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]][1] >= 1){
      poly@polygons[[i]]@Polygons  <- poly@polygons[[i]]@Polygons[bigpolys[[i]]]
      poly@polygons[[i]]@plotOrder <- 1:(length(poly@polygons[[i]]@Polygons))
    }
  }
  return(poly)
}

reduced.shapefile <- getSmallPolys( simplified.shapefile , 0.01 )

This worked perfectly as long as a stayed with the basic R plot() function.

However, shifting to the Leaflet package caused problems. While the other SpatialPolygonsDataframes would work fine (albeit very slowly) in Leaflet, the reduced.shapefile would not.

leaflet(data = reduced.shapefile) %>%
  addPolygons(fillColor = "blue",
              fillOpacity = 1,
              color = "transparent",
              weight = 1)

Running that code would give the error: Error in pgons@Polygons[[index]] : subscript out of bounds.

Can anyone see the specific cause of the problem, and a way that I can restructure or adjust the function (getSmallPolys) such that the new SpatialPolygonsDataframe will work in Leaflet and not just plot()?

I've included some additional notes below:

If anyone's interested, the SA2 ESRI shapefile I'm using can be downloaded here: http://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/1270.0.55.001July%202011?OpenDocument
Examples of the specific subgroups of polygons that cause the problems include reduced.shapefile[156,] and reduced.shapefile[981,]
I suspect that the problem relates to a part of Leaflet trying to reference an element of i in  reduced.shapefile[981,]@polygons[[1]]@Polygons[[i]] that is larger than the number of polygons in the reduced shapefile.


texts: Given a list of surveyor points, what is the algorithm used to simplify the points and retain the terrain characteristics?

texts:
Given a terrain and surveyor points sampled from from it or from its contours, is there any algorithm that one can use to simplify the points ( i.e., reduce the surveyor points number) and at the same time retaining the terrain characteristics?

The reason I ask this is because the surveyors may take too much redundant points at flat terrain, and these points serve nothing but to slow down my volume calculation process, so I want them to get filtered out.


texts: To load a raster layer in ArcMap with the C# application

texts:
I have a Python script which import a layer in ArcMap, when I execute it in Python console of ArcMap  but now I have to perform the same work with  C# program.
Python script for importing a raster layer is as:

import arcpy
from arcpy import *
arcpy.env.workspace ="D:/Data/Old_Data"
raster_filename ="output.tif"
arcpy.MakeRasterLayer_management(raster_filename, "tmpLyr")

Can creating a toolbox for the script help for this?


texts: Listen to file dropped in openlayers 3 Drag and Drop interaction

texts:
I'm trying to display an SHP file into openlayers 3 map. I found a library to read both .shp and .dbf files, so I can read features and attributes from the files. Now I want the client to drop those files to the map using DragAndDrop interaction.

Is possible to do this?

Is there an event similar to bellow example?

dragAndDrop.on('file-drop', function(event){
  // Do something with event.file or event.someKey
});

In this event (if exist) should I parse the files, get features and add to map manually?

Also if this kind of event exist what would happen with addfeatures event.

Will be triggered??


texts: Why is 'Create Geospatial PDF.' greyed out?

texts: I am using QGIS to create a map that can be used in avenza. When I try to export as a PDF it does not allow me to tick the box of Create geospatial PDF. I am using Mac Catalina and the current version of QGIS is 3.14.


texts: arcpy.MakeTableView returns empty result object on second iteration of loop

texts:
I've encountered a very strange problem involving a geoprocessing tool in ArcPy that I have spent too much time trying to solve on my own.
Specifically, I have a process that iterates through a list of folders with images, adds those images to an existing mosaic dataset (in a file geodatabase), and then updates metadata for the added raster within the mosaic dataset attribute table, based on metadata stored in an associated JSON file.
Earlier in this workflow, a previously-run tool generates folder names (based on an image name obtained from the JSON file), and if the process finds a folder name will be duplicated, the folder name is appended with a numerical suffix (i.e. '..\20200509_135639_1010', '..\20200509_135639_1010_01' ).
All raster files are added to the mosaic dataset w/o issue. However, when I try to update the metadata associated with each rasterin the mosaic dataset attribute table (again based on values stored in an associated JSON file) the problem occurs.
I attempt to create a temporary Table View for one raster record, based on a query in which I select the newly-added raster based on the image's "product name" in the mosaic dataset attribute table. This works for the first image in a list, but NOT for the subsequent images which have the appended suffix (i.e. "_01", "_02", etc.)
No matter how I change the values associated with the suffix ("01" or "_1", etc.) the arcpy.MakeTableView_management always return an empty result object on all subsequent "suffixed" raster names.
Here is the offending code snippet that returns an empty result on the "suffixed" raster names:
Example variable values

srcMosaic = '\\\\SERVER_NAME\\DIR\\SUBDIR\\SourceMosaics.gdb\\PlanetScope'
whereClause = "ProductName = '20200509_135639_1020_01'

def updateMetadataFromUser(srcMosaic, whereClause, collectionName, collectionType, imageType, source, sensor, year, startDate, endDate, resolution_meters, loadedBy):
    metaTable = arcpy.MakeTableView_management(in_table=srcMosaic, out_view="tmpMetaTable", where_clause=whereClause)

    #Populate the table with attributes collected from getRasterCollectionInfo and getRasterInfo
    metaCount = arcpy.GetCount_management(metaTable)[0]

The metaCount = arcpy.GetCount_management(metaTable)[0] returns an error when the metaTable result object is empty.
Here's an example of the image names:

'20200509_135639_1020' (metadata fields updated successfully in mosaic attribute table)
'20200509_135639_1020_01' (returns empty result from MakeTableView)

Another bizarre aspect to this is that if I just run the tool using only the image with the suffix (for instance 20200509_135639_1020_01), then arcpy.MakeTableView find the record and returns a result.
Any ideas?


texts: Seeking QGIS 3.0 ways to create browser based 3D scenes like qgis2threejs plugin does?

texts:
I want to create a browser-based 3D scenes and use Github to share it. However, the qgis2threejs plugin is not available in QGIS 3.

Is there any alternative way to create similar browser-based 3D scenes?


texts:
I am making an interactive webmap with QGIS2WEB. Normally when I create layers in QGIS then open QGIS2WEB the layers open up in the preview window, then appear in the website.

However, I now have some custom layers and not only will they not appear in the preview window or the created web map, if I include them at all the whole map shows up blank.

I tried this with layers whose CRS was EPSG 4326 as well as 3857, neither worked.

I tried opening up some random shapefiles downloaded from the web and they worked fine, but once I added my layers, any of them, the whole map turned up blank.

I know there are some limitations to what QGIS2WEB can display, but these particular .shp layers are fairly basic and work fine elsewhere...

Any advice?


texts: qgis-earthengine: using a local qgis layer feature (e.g. shapefile polygon) as a roi

texts:
I wanted to use a local shapefile polygon as roi; i.e.

roi = ee.FeatureCollection("path/to/the/qgis_project_name/polygon_name")

but, I got:

  Collection.loadTable: Collection asset
  'path/to/the/qgis_project_name/polygon_name' not found.

I tried to export a shapefile in Google drive but I also got the same error.


texts: Selecting the active layer in Layers tree in QGIS

texts:
How do I programmatically (in QGIS with Python) set the position of the active layer and show it in the legend tree.
Specifically if the layer is in a collapsed group or is it far below the list?


texts:
First of all, I have set the WapMapServiceImageryProvider with Geoserver deploy radar image to the Cesium, The result like this:

The code is :

function addRadarHKZindex(physicalHeight){
  if(!activeRadarHeight.hasOwnProperty(physicalHeight)){
    var geoWmsLayer = new Cesium.WebMapServiceImageryProvider({
      url : geoserver_url,
      layers: 'layer_name',
      parameters : {
            service     : 'WMS',
            version     : '1.1.0',
            request     : 'GetMap',
            styles      : radarHeight[physicalHeight],
/* The Geoserver determines style to define the radar height */
            format      : 'image/png',
            transparent : true,
            srs         : 'EPSG:4326',
            width       : '545',
            height      : '768'
        }
    });
    activeRadarHeight[physicalHeight] = imageryLayers.addImageryProvider(geoWmsLayer);
    scene.requestRender();
  }
}

But the radar image is on the ground with 0 ft, I would like to display WebMapServiceImageryProvider with Physical Height on Cesium.

The idea just like this

I have tried use Rectangle + height with WebMapServiceImageryProvider but not working, here's the code

function addRadarHKZindex(physicalHeight){
  if(!activeRadarHeight.hasOwnProperty(physicalHeight)){
    var tmp_rectangle = Cesium.Rectangle.fromDegrees(111.32, 19.8, 117.24, 25.9);
    tmp_rectangle.height = physicalHeight;

var geoWmsLayer = new Cesium.WebMapServiceImageryProvider({
  url : geoserver_url,
  layers: 'layer_name',
  rectangle : tmp_rectangle, // Add physicalHeight into rectangle
  parameters : {
        service     : 'WMS',
        version     : '1.1.0',
        request     : 'GetMap',
        styles      : radarHeight[physicalHeight],

/* The Geoserver determines style to define the radar height */
            format      : 'image/png',
            transparent : true,
            srs         : 'EPSG:4326',
            width       : '545',
            height      : '768'
        }
    });
    activeRadarHeight[physicalHeight] = imageryLayers.addImageryProvider(geoWmsLayer);
    scene.requestRender();
  }
}

However, the result is not working, the result is showing grey rectangle:

so How to display WebMapServiceImageryProvider with Height on Cesium correctly?


texts: Comparing two kernel density rasters for areas of interest?

texts:
The background to my question is that I have unit sales for Product A aggregated by city. I also have unit sales for Category A aggregated by city.

The goal is to compare areas where Product A performs better than expected given performance of Category A. The reason for comparison is that unit counts tend to follow general population distributions (ie. at face value both data sets show same big cities are the primary markets).

What I did was take each data set and normalize it on z-score. I then computed two kernel density rasters, one for each set. Now my question is if there is a best way to compare the two given my objective?

Can I just use a raster calculation [Product A] - [Category A] to create a raster to show areas where the product outperforms the category? Is that a valid approach or is there a better way?

Also I noticed that the min/max value range computed for each kernel density is different for the two data sets. Does this matter for comparison?

Product A: -0.0055 to +0.0253
Category A: -0.0087 to +0.0683

Update
The more I play with this, I'm not getting the results I need. I can't use division operators in the raster calculation because near-zero demonstrators are astronomically inflating the numbers. I believe I would have to re-normalize each of the two kernel densities first?


texts: ArcGIS Viewer for Flex supported EPSG (wkid)

texts:
I am trying to "force" ArcGIS Viewer for Flex 2.1 to retrieve WMS layer using an EPSG, i.e. 3308.

I tried setting

 <map wkid="3308">

However the ArcGIS Viewer for flex always defaults to 4326, unless I include the default basemap, in which case the SRS in WMS http request is 12100 (or something similar)

<basemaps>
        <layer label="Streets" type="tiled" visible="false"
               url="http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer"/>

So how do I display other WMS with EPSG 3308 ?


texts: Deleting Feature Class from Geodatabase using list in ArcPy?

texts:
I am trying to delete several feature classes from a file geodatabase using values stored in a user defined list.

The script runs without error but it does not delete the feature classes.

Any suggestions to what I am doing wrong (using ArcMap 10.3)?

import arcpy

arcpy.env.workspace = r"C:\Data\Temp.gdb"

fc_Delete = ["fcOut1a","fc_Out2a","fc_Out3a"]

for fc in fc_Delete:
    if arcpy.Exists(fc):
        arcpy.Delete_management(fc,"")


texts: Using PrintTask a Map with Graphics symbols not solid

texts:
I have a Web map using JavaScript API.  One of the layers is a Graphics Layer that I add graphics to.

The line symbol style for these graphics I add is STYLE_SHORTDOT.  When I have a graphics drawn on the map, the PrintTask task refuses to export the map.  (It works as long as that type of graphic is not on the map.)   All I have to do to get the PrintTask to like my graphics is use STYLE_SOLID.

Has anyone else seen this behavior?  Am I condemned to use only solid line graphics if I want to export my map?


texts: How to programmatically Summarize an ITable?

texts:
Is there a way in ArcMap plugin to make Summarize-function like functionality programmatically?

There is a feature layer, lets say Building, And it has fields like OID, FloorNumber, Color, Height etc... So I want summary like:

Floors TotalHeight Heightest
3      120         20 ( 4 floor buildings in total have 120 meters height )
4      ...         ...


texts:
I've got a bunch of data from LINZ which includes rocky areas. Is it possible to style them as they do (see picture)? It's a fairly standard style in large scale maps so I hoped it would be simple!

I'm using QGIS but I'm relatively new to it


texts:
I am new to QGIS (2.6.1) and am having trouble adding a new linear feature to an existing layer.

I have set tolerances and digitized the vertices of the line and it terminates with a large green X at the terminal node.  But any attempt to select the feature results in 0 features found. And if I turn off the layer, this solitary line remains visible.  What layer is the line associated with?  What am I missing?


texts: How do I upload a Featureclass to a Geodatabase in C#?

texts:
I have a featureclass with Arcobjects that came from a Shapefile and I wish to import it into a geodatabase. I have the geodatabase connection, I just don't know how to load the feature class to a Geodatabase.  Any help would be amazing!

I'm looking for a code example seeing I can't find it by googling. Hopefully in C#, but I don't mind reading through VB if needed.

EDIT:

I wanted to add the Answer thanks to George so here is a code snippet with comments.

 ESRI.ArcGIS.Geodatabase.IWorkspace pWorkspace;
        ESRI.ArcGIS.Geodatabase.IWorkspaceFactory pWorkspaceFactory;
        ESRI.ArcGIS.Geodatabase.IFeatureWorkspace SdeFeatureWorkspace;

        pWorkspaceFactory = new FileGDBWorkspaceFactory();
        pWorkspace = pWorkspaceFactory.OpenFromFile(@"C:\tests\test.gdb", 0); //Opens the Geodeabase and creates a connection.
        SdeFeatureWorkspace = (ESRI.ArcGIS.Geodatabase.IFeatureWorkspace)pWorkspace;

        IWorkspaceFactory workspaceFactoryShape = new ESRI.ArcGIS.DataSourcesFile.ShapefileWorkspaceFactoryClass(); //creates a shape factory to pull the shape and create the feature class.
        IFeatureWorkspace featureWorkspace = (IFeatureWorkspace)workspaceFactoryShape.OpenFromFile(@"C:\Users\Desktop\Newfolder", 0); //Opens the Shape file.
        IFeatureClass featureClass = featureWorkspace.OpenFeatureClass("al012009.001_5day_lin"); //opens the featurclass from the shapefile

        SdeFeatureWorkspace.CreateFeatureClass("al012009001_5day_linnn", featureClass.Fields, featureClass.CLSID, featureClass.EXTCLSID, featureClass.FeatureType, featureClass.ShapeFieldName, null); //creates a feature class in the geodb.

        Geoprocessor gp = new Geoprocessor(); //starts a GEOprocess.
        gp.SetEnvironmentValue("workspace", @"C:\Tests\test.gdb"); //sets the geodb to use

        ESRI.ArcGIS.DataManagementTools.Append  append = new ESRI.ArcGIS.DataManagementTools.Append(); //this tool uploads the feature class
        append.inputs = featureClass;
        append.target = @"C:\tests\test.gdb\al012009001_5day_linnn";

        gp.Execute(append, null);


texts:
I have downloaded a map file which I'm using as a texture on a plane. The file is a proper map file with projection in a prj-file and a tfw-file. The prj file contains:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]]

And the tfw-file contains:

      0.00000000000000
      0.00000000000000
      -0.03333333333333
      -179.98333333333333
      89.98333333333333

The original file is 10800x5400 pixels, and I've scaled it down to 2048x1024 pixels. So, how do I calculate the long/lat pair if someone clicks on point [712, 133] on my scaled down image?

Bonus things I'm wondering about:

The image will be rescaled several times, but will keep the aspect ratio, is that a problem?
How do I do this in reverse? IE if I want to find the x,y
pair from a long/lat pair.

I've spent hours trying to read in random places about GIS, maps and whatnot but it seems to be a very steep learning curve and I'd be happy to take a shortcut as in the whole it's a fairly minor part of what I'm trying to do. Pseudo code, c#, some other language doesn't really matter, but I'm surprised how hard it has proven (for me at least) to find such hands on advice.

Sorry about the cross post, asked here originally but based on the comment I figured it might make more sense to ask here


texts: Determining vegetation health on temporal scale using mean NDVI?

texts:
I’m calculating NDVI for Bryce Canyon National Park over many years.  I want to find if the vegetation in this area is becoming less healthy due to drought and heat stress over the years. When I calculate mean and max values, the mean value seems to be far off.  For instance, for one of the years: Max NDVI: 0.70, Mean NDVI: 0.05.

I'm consistently getting numbers in this low range (0.03 - 0.05).
I’m thinking it has to do with the majority of the map being soil.

So, I'm wondering if using the mean NDVI is a good approach for seeing if the health of vegetation in this area is declining or getting better?

I've tried using SAVI, but this doesn't change the mean and max by much.

Also, I'm aware the NDVI anomaly approach, but I'm unsure as to how to calculate this or if this would be a better approach.


texts: Discrepancy between shape area and data from municipality

texts:
I am attempting to map several parcels just to show very basic dimensions (frontage, etc.) and parcel area. I have the parcel area from the municipal assessor's office in acres. And I have a parcel shapefile from my state's GIS entity (MassGIS) that I believe is drawn from the municipal data. There is a "Shape Area" field in the parcel shapefile attribute table, which appears to be in sq meters. I added a field for Area in sq feet, and did "calculate geometry". My problem is that the values in the attribute table for the parcel shapefile are significantly different from those provided by the assessor's office, and I am not sure whose numbers to trust. I need the dimensions of parcels to be correct, and I'm worried that the data in my shapefile is either incorrect or being skewed somehow, which will result incorrect dimensions.

Has anyone else run into this problem? Does this sound like a projection issue, or an issue with the data itself?


texts:
I have data captured with a GPS tracking unit. Some track have bad/cold starts and the only way to determine the track quality is by using the altitude. Unfortunately I don't have the raw GPS data and only an array or Lat, Long & Altitude.

A typical bad start will have an altitude that start at negative with a steep ramp up to the proper altitude.

Is there any maths gurus and/or libraries (code) out there that can help me with an algorithm to detect and clean my track data?


texts: I'd like to learn how to use a neural network in GIS. I surfed in the internet to find good resources or a tutorial about this but I couldn't find any helpful resources. Are there any good references or tutorials about using neural networks in GIS? I'd like to use an artificial neural network (ANN) in mapping hazard zones, for example flood zones or landslides. I use ArcGIS 10 to process the layers but don't know how can  use an ANN.


texts: Can't add map service (map cache) to mxd

texts:
When trying to add a map service (in this case a map cache service) to an ArcMap project, I keep getting an error that says:

Could not add specified data object to the map.
Could not connect to the Map Server.

Even though in Catalog I can connect to the service perfectly fine and open it as an mxd project, I can't add it to another project. The goal is to use the map cache as a basemap for editing natural features from an aerial photo. Any help would be greatly appreciated.

Using ArcMap 10, ArcServer 10

Thanks!


texts: Coordinates translation issue (from EPSG:23030 to EPSG:4326)

texts:
I'm trying to figure out how to make some translations from UTM/ED50 to WGS84.

The reason I need to do them is that my data source's coordinates are in EPSG:23030 format, and I need to represent them in a map, currently using leaflet, which needs every coordinates to be in WGS84.

I've been using Postgis' raw sql queries to transform coordinates, as follows:

SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT(581087 4796507)', 23030), 4326));

That query returns me the following:

"POINT(-2.0013194711113 43.3151520669395)"

As I'm developing a web project using html+js and so on, I needed some library to make translations. So I found proj4js, and I can make translations as follows:

proj4.defs([
    [
        'EPSG:4326',
        '+title=WGS 84 (long/lat) +proj=longlat +ellps=WGS84 +datum=WGS84 +units=degrees'
    ],
    [
        'EPSG:23030',
        '+title=ED50 (UTM) +proj=utm +zone=30 +ellps=intl +units=m +no_defs'
    ]
]);

var coord = proj4('EPSG:23030', 'EPSG:4326').forward([581087,4796507]);

holding coord the following coordinate:

43.315428738109816, -2.000074748128657

So the question is why am I getting different coordinates, knowing that I'm -seemingly- using the same coordinate systems to make the translation in PostGis and in Proj4js? Is there something wrong?

As far as I know, the right conversion corresponds to PostGis.


texts: Why is Clip (coverage) in ArcMap not working?

texts:
Clip Analysis and Clip Data Management work fine, but I cannot open the tool or the help menu for Clip Coverage.

Does this sometimes happen because the tool requires an extension to run?


texts: Labeling Many-to-many relationships ArcMap 10

texts:
please i have Many-to-many relationships in arcMap like :

and i want to show labels in parcel  using ArcMap 10 :

field to show is : Zoning and Name
 -

do you have any idea how to do it ?

Edit
 -


texts: Choropleth Map interactively changing style

texts:
I'm trying to create a interactive choropleth web map with qgis2web, but having some difficulties. I know how to create a web map and make some information pop up as you click an especific feature.

But I also want to change the Style of the whole map as you click or put the mouse on a certain feature. The real change would be on the column selected.
An example:

If you put your mouse on "Bangu" feature, the "Bangu" column would be selected and the map shown would be like this:

But if you put your mouse on another feature, like "Flamengo", the "Flamengo" column would be the selected one and the map shown would be like this:

I've thought about using leaflet, but not really sure of how.

Any ideas, suggestions, orientations, advices, informations?


texts: Altering existing polyline shapefile that displays contour lines using ArcGIS Desktop?

texts:
I have a polyline shapefile that displays contour lines. It has contour lines that have a interval of 10 m.

Is there any way to change the interval to 500 m in ArcGIS?


texts:
I need to suggest a tool to build a map together, ie. let multiple users add POIs and routes on a map.

I know about OSM-based Umap but are there other online tools I should know about?


texts: How do I create a common symbology scale for multiple rasters in QGIS/ArcGIS?

texts: I have twelve global surface temperature rasters and each has a slightly different temperature range depending on the time of year. I would like to create a common scale i.e. -80c to 80c and so all the colours in each raster match each other. I am happy it do this in either QGIS or ArcGIS.


texts: Adding drop down option for layers in Python script tool?

texts:
I've code that will add fields to feature classes.  I've parameterized the script so the user can navigate to the location of the feature classes to be edited.

If the feature class is already open and in my data frame I want the option for the drop down to pick feature classes that are already open rather than having to navigate to the location of the files.

Does that require code being added to my script or is that something I can do in the script properties?

My code looks like this currently.

import arcpy

fc = arcpy.GetParametersAsText(0)

# Process: Add Field
arcpy.AddField_management(fc, "DATE_MOD", "DATE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: Add Field
arcpy.AddField_management(fc, "DATE_ADDED", "DATE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: Add Field
arcpy.AddField_management(fc, "MOD_BY", "TEXT", "", "", "50", "", "NULLABLE", "NON_REQUIRED", "")

# Process: Add Field
arcpy.AddField_management(fc, "ADDED_BY", "TEXT", "", "", "50", "", "NULLABLE", "NON_REQUIRED", "")


texts: Is there a simple interface from which I can extract average summer temperature and humidity for a list of latitude and longitudes?

texts:
I am interested in finding the summer average (either climatalogical mean or for a specific year) relative humidity and temperature for a set of globally distributed latitude and longitude coordinates.

What I need are mean values for June, July, and August.

I would prefer to use actual observed values, but a reanalysis (such as NARR or CRU-NCEP would suffice). Is there a web-based interface from which I can obtain this information?

The closest that I have come to this

Wolfram Alpha, but they only provide rainfall and temperature, not relative humidity
worldclim.org looks promising, but requires me to download all of the data before using another program to extract information; it also does not have relative humidity.


texts: Enabling ability to change basemap in Carto?

texts: I have a carto web map that I would like to be able to toggle between the Carto basemap and satellite imagery. Similar to how this is easily done in ArcGIS Online with "Basemap Gallery" widget.


texts: Length of intersections from a linestring and a grid shapefile by using Python GeoPandas & QGIS

texts:
I am using Python GeoPandas and I have two shapefiles. One of them is a grid file containing equal polygons. The other shapefile is a linestring file with several lines and its Attributes.

The idea is to do a multiplication with the intersection length of each individual line within a grid cell and finally sum up all multiplication results which belong to this certain grid cell. Each row is one grid cell in the grid shapefile. I need the length of each intersection which is inside a grid cell. So there might be only one grid cell, but multiple intersections within this certain grid cell.

This simple multiplication should be like this:

df['new_column_result']= line['intersection_length_inside _cell_i']*line['factor_x']

And afterwards use this multiplication to calculate the sum of all df['new_column_result'] that are in the cell.

Additionally, the grid cells which are not touching any lines at all should be removed and the two shapefiles should be merged as an output. This part I managed to do by using the sjoin function of GeoPandas, but I am not sure if it is suitable for my goals. Sadly at the moment I get the same grid cell (with the same grid cell ID) multiple times. I guess it is produced for each "intersection_length_inside _cell_i".

For joining this two dataframes(or shapefiles) I use:

intersections = geopandas.sjoin(gridshapefile, lineshapefile, how="inner", op='intersects')


texts: Why does reprojection with gdalwarp change pixel size?

texts:
Is it okay that reprojection from UTM (EPSG:32612) to Web Mercator (EPSG:3857) with gdalwarp changes pixel size? Why does this happen? Does that mean that unit is not meters? If so, than how to convert pixel size from UTM units to Web Mercator units?

Command:

gdalwarp -t_srs "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs" LC08_L1TP_047011_20170518_20170519_01_RT_B4.TIF output.tif

Input image is Landsat 8 from AWS:
http://landsat-pds.s3.amazonaws.com/c1/L8/047/011/LC08_L1TP_047011_20170518_20170519_01_RT/LC08_L1TP_047011_20170518_20170519_01_RT_B4.TIF

Original image info:

$ gdalinfo LC08_L1TP_047011_20170518_20170519_01_RT_B4.TIF
Driver: GTiff/GeoTIFF
Files: /home/user/Desktop/Landsat_sample/LC08_L1TP_047011_20170518_20170519_01_RT_B4.TIF
Size is 8451, 8491
Coordinate System is:
PROJCS["WGS 84 / UTM zone 12N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-111],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32612"]]
Origin = (488985.000000000000000,7851615.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  488985.000, 7851615.000) (111d17'58.67"W, 70d46' 6.92"N)
Lower Left  (  488985.000, 7596885.000) (111d16' 8.84"W, 68d29' 3.48"N)
Upper Right (  742515.000, 7851615.000) (104d25'49.46"W, 70d39' 3.78"N)
Lower Right (  742515.000, 7596885.000) (105d 5'39.63"W, 68d22'48.67"N)
Center      (  615750.000, 7724250.000) (108d 1'24.05"W, 69d36' 5.29"N)
Band 1 Block=512x512 Type=UInt16, ColorInterp=Gray

Reprojected image info:

$ gdalinfo output.tif
Driver: GTiff/GeoTIFF
Files: out.tif
       out.tif.aux.xml
Size is 8902, 8880
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-12389818.176861140877008,11323689.164804801344872)
Pixel Size = (85.895238477817642,-85.895238477817642)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-12389818.177,11323689.165) (111d17'58.67"W, 70d46' 7.79"N)
Lower Left  (-12389818.177,10560939.447) (111d17'58.67"W, 68d22'48.98"N)
Upper Right (-11625178.764,11323689.165) (104d25'50.73"W, 70d46' 7.79"N)
Lower Right (-11625178.764,10560939.447) (104d25'50.73"W, 68d22'48.98"N)
Center      (-12007498.470,10942314.306) (107d51'54.70"W, 69d36'28.82"N)
Band 1 Block=8902x1 Type=UInt16, ColorInterp=Gray
  Min=0.000 Max=52607.000
  Minimum=0.000, Maximum=52607.000, Mean=17537.063, StdDev=16780.853
  Metadata:
    STATISTICS_MAXIMUM=52607
    STATISTICS_MEAN=17537.062597405
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=16780.852772678


texts: Aggregating multiple layers while preserving values using ArcMap?

texts:
I am using ArcMap 10.5 to display the distribution of Revenues by Region, for each product I have defined.

Each product has a layer of its own consisting of a collection of polygons, and each polygon contains attribute fields for revenues over a period of 10 years (e.g. Revenues2005, Revenues2006, etc)

I am trying to create a kind of cumulative heatmap.

For example, let's take two polygons with the following fields, each in a separate layer, but that overlap exactly:

Layer1

Rev1 | Rev2 | Rev3

5000   | 6000 | 7000

And

Layer2

Rev1 | Rev2 | Rev3

1000   | 2000 | 2500

Would it be possible to combine the above layers, such that when the polygons overlap, I would obtain a new layer with a polygon of the following attributes:

Layer3

Rev1 | Rev2 | Rev3

6000   | 8000 | 9500

So far the (tedious) solution I found is to combine the fields pairwise and use a python script to add the revenues for each year one by one.


texts: How to get edge length(distance) from pgrouting DB

texts:
I have implemented pgRouting workshop and it's working fine. I need to get every edges length in kilometers to do some calculations I studied the ways table. It contains a length attribute, but it seems to be not what I exactly need.

Is there any place or way to get the length of each and every edge in kilometers or meters.


texts: Using Leaflet Search plugin with GeoJSON layer (L.Control.Search)

texts:
I have a map with a single GeoJSON layer that I referenced in my script (using the AJAX plugin). Now I am trying to add search functionality to the map so that the user can search for features by name. The plugin I am using is the Leaflet Search plugin.

However, I cannot get the search to perform, as it's giving me the "Location not found" error, whenever I type text into the search box. Here is my code:

searchlayer = L.geoJSON.ajax('data.json').addTo(map);
map.addControl(new L.Control.Search({layer:searchlayer});

What am I doing wrong here? I assume that I'm not properly linking the layer and the field to be searched to the search control, but how do I do that?

Found that the issue was not adding the option for the control to search the geojson attribute. The default column that the search button works on is "title". I added the option of propertyName: 'Location' to L.Control.Search and it works perfect now!


texts: Images do not appear on the map in Google Earth Engine, no bands to visualize

texts:
I have a code that shows how to make a composite image for every quarter from 2013 to 2018. However when I tried to print some of them, I couldn't see any images or layers in the map below.

The code which I got from https://code.earthengine.google.com/f7a4cf857dca37225213c43a879b3f15
// get start and end date and define the number of quarterly composites
var start = ee.Date('2013-01-01');
var end = ee.Date('2019-01-01');
var numbQuarters = end.difference(start, 'month').divide(3).ceil();

// make a composite image for every quarter
var quarterlyImages = ee.ImageCollection.fromImages(
            ee.List.sequence(1, numbQuarters).map(function(quarter){
              var startTemp = start.advance(ee.Number(quarter).subtract(1).multiply(3), 'month');
              var endTemp = start.advance(ee.Number(quarter).multiply(3), 'month');
              var image = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
                                .filterDate(startTemp, endTemp)
                                .filterBounds(geometry)
                                // define a way to composite an image
                                .median(); // .mean(); .mossaic(); .max(); etc

              return image.set('system:time_start', startTemp.millis(),
                               'system:time_end', endTemp.millis());
            }));

// print some stuff
Map.addLayer(quarterlyImages)
print('startDates', quarterlyImages.aggregate_array('system:time_start'))
print('sendDates', quarterlyImages.aggregate_array('system:time_end'))


texts:
I'm trying to run a Python program which needs to import

from osgeo import gdal, ogr
import osgeo.ogr, osgeo.osr
from osgeo.gdalconst import *

But when I run it outputs

Traceback (most recent call last):
  File "./lisbon_P2M.py", line 8, in <module>
    from osgeo import gdal, ogr
ImportError: No module named 'osgeo'

Meanwhile, I tried to install osgeo but it pops this error:

sudo apt-get install libgeos-dev

Reading package lists... Done
Building dependency tree
Reading state information... Done
Some packages could not be installed. This may mean that you have
requested an impossible situation or if you are using the unstable
distribution that some required packages have not yet been created
or been moved out of Incoming.
The following information may help to resolve the situation:
The following packages have unmet dependencies:
 libgeos-dev : Depends: libgeos-c1v5 (= 3.5.0-1ubuntu2) but 3.5.1-3~xenial0 is to be installed
E: Unable to correct problems, you have held broken packages.

Any hints on how this can be solved ?


texts: Web Appbuilder Widget refreshes and changes popup visibility to hidden

texts:
I am trying to update an app used from Web Appbuilder 1.2 to WAB 2.11.  I cannot get the nexrad widget/map to behave the same in the newer version of WAB.  In 1.2, the widget refreshes (1,2,5,or 10 second intervals).  In 2.11, every time the Widget refreshes, my layer popups (when I click a feature on the map) disappear.  Using developer tools, the styles read as follows:

.esriPopupHidden,
.esriPopupHidden * {
    visibility: hidden;
    opacity:  0;

I can uncheck these with the developer tools and the popup is visible again.  I tried adding

.esriPopupHidden,
.esriPopupHidden * {
    visibility: visible !important;
    opacity:  1 !important;

to the styles.css file but that doesn't seem to work.  Is there somewhere in the CSS or JS where I need to make a change to decouple these actions (Nexrad refreshes and layer popups)?  There is no obvious location with the div id or an esri.css file.


texts:
Is there any way to append a shapefile in a one-to-many relationship with its attributes that R can read? I understand the fix that resides in ArcInfo by creating through MakeQueryTable, but as this information is being analyzed in R -- that is unneccessary.

Really, I am trying to understand R's capabilities with one-to-many relationships as well as the limitations it has utilizing shapefiles.


texts: esri.tasks.GenerateRendererTask for non 10.1 ArcServer

texts:
I wish to use the same features as the esri.tasks.GenerateRendererTask but it seems that it can be only used on ArcServer 10.1. Is it possible to achieve the same dynamic natural jenks breaks in the ArcServer 10?

So I will be having a ArcGISDynamicMapServiceLayer and I have to generate the renderer on the fly and the values are also dynamic and not in a fixed range. esri.tasks.GenerateRendererTask takes care of it beautifully, is there any other solution to this problem?


texts: Export data from time-series image collection to a table in Google Earth Engine

texts:
On Google Earth Engine, I'd like to accomplish two tasks:

extract data from every pixel, across a time-series, within the given geometry bounds and save it as a table.
I need a column which prints the date of the image.

For the first part here is my code snippet (borrowed from : Extract complete pixel values inside a geometry) ---

// the geometry bounds
var rectangle = ee.Geometry.Rectangle(96.01669, 18.52621, 96.04819, 18.49634);

//the time series dataset
var dataset = ee.ImageCollection("MODIS/006/MOD13Q1")
.filterDate('2000-01-01', '2000-03-31')
.filterBounds(rectangle)
.select('NDVI', 'SummaryQA');

// Sort the dataset
var sorteddataset = dataset.sort('system:time_start', true);

var first = sorteddataset.first();

// generate a new image containing lat/lon of the pixel and reproject
// it to MODIS projection
var coordsImage =
ee.Image.pixelLonLat().reproject(first.projection());

var joinedImage = coordsImage.addBands(first);

var valuesList = joinedImage.reduceRegion({
reducer: ee.Reducer.toList(4),
geometry: rectangle
}).values().get(0);

valuesList = ee.List(valuesList) // Cast valuesList

var myFeatures = ee.FeatureCollection(valuesList.map(function(el){
el = ee.List(el) // cast every element of the list
var geom = ee.Geometry.Point([ee.Number(el.get(0)),
ee.Number(el.get(1))])
return ee.Feature(geom, {'NDVI':ee.Number(el.get(2)),
'SummaryQA':ee.Number(el.get(3))})

}))

Export.table.toDrive(results,
"data",
"ndvi-testing",
"NDVI-SummaryQA",
"CSV");

How do I modify the code to include the entire time-series instead of just the first image?
And how can I include a command that prints the image date?
Below is the image of the table I get from the given code.


texts: How to calculate straight line distance from spherical distance

texts: I have longtitude and lattitude coordinates of one point and spherical distance between this and another point on Earth in kilometers (the shortest distance over the earth’s surface). I need to calculate what is the distance between these points in straight line. Can anybody help?


texts: Display the output generated by saga ordinary kriging on to QGIS from Python Plugin

texts:
I wanted to call kriging algorithm from a python plugin and displays it's output on QGIS desktop.

I have successfully generated the output of the kriging algorithm by unable to find the steps which will help to display the output on QGIS desktop as far as I know there is no help in document(I would be happy if anyone can find and suggest me.)

Can anyone suggest me what should I do to solve this problem.
Below is the lines which I have used to generate the output.
OS:UBUNTU

import processing
Output = processing.runalg('saga:ordinarykriging',
                           "/home/harsh/Documents/IIC/Krigging_test/demo.shp", #POINTS <ParameterVector>
                           "cost", #FIELD <parameters from POINTS>
                           0, #TQUALITY <ParameterSelection>, [0] standard deviation
                           True, #LOG <ParameterBoolean>
                           True, #BLOCK <ParameterBoolean>
                           100,  #DBLOCK <ParameterNumber>
                           -1,   #VAR_MAXDIST <ParameterNumber>
                           100,  #VAR_NCLASSES <ParameterNumber>
                           1,    #VAR_NSKIP <ParameterNumber>
                           "a + b * x", #VAR_MODEL <ParameterString>
                           "0, 1, 0, 1", #OUTPUT_EXTENT <ParameterExtent>
                           100.0, #TARGET_USER_SIZE <ParameterNumber>
                           0, #TARGET_USER_FITS <ParameterSelection>, [0] nodes
                           0, #SEARCH_RANGE <ParameterSelection>, [0] local
                           1000, #SEARCH_RADIUS <ParameterNumber>
                           0, #SEARCH_POINTS_ALL <ParameterSelection>,  [0] maximum number of nearest points
                           4, #SEARCH_POINTS_MIN <ParameterNumber>
                           20, #SEARCH_POINTS_MAX <ParameterNumber>
                           0, #SEARCH_DIRECTION <ParameterSelection>, [0] all directions
                           None, #PREDICTION <OutputRaster>
                           None) #VARIANCE <OutputRaster>

print Output
{'VARIANCE': u'/tmp/processing5f15099df76849648e0f59eea0f3d625/00565712f1274a6082cad84e1a5a944c/VARIANCE.sdat', 'PREDICTION': u'/tmp/processing5f15099df76849648e0f59eea0f3d625/32a911a50c3e480cbf9418a8593d38ca/PREDICTION.sdat'}

[Solved]

ext = QgsVectorLayer( '/home/harsh/Documents/IIC/Krigging_test/demo.shp', '', 'ogr' ).extent()
        xmin = ext.xMinimum()
        xmax = ext.xMaximum()
        ymin = ext.yMinimum()
        ymax = ext.yMaximum()
        coords = "%f,%f,%f,%f" %(xmin, xmax, ymin, ymax) # this is a string that stores the coordinates

        Output = processing.runalg('saga:ordinarykriging', "/home/harsh/Documents/IIC/Krigging_test/demo.shp", "Z", 0, True, True, 100, -1, 100, 1, "a + b * x", coords,100.0,0,0,1000, 0, 4,20,0, None, None)
        print Output
        prediction = QgsRasterLayer(Output['PREDICTION'],
                                    'prediction')

        QgsMapLayerRegistry.instance().addMapLayer(prediction)

        variance = QgsRasterLayer(Output['VARIANCE'],
                      'variance')

        QgsMapLayerRegistry.instance().addMapLayer(variance)


texts:
I'm using ArcGIS 10.1 with ArcPad 10.0.  I am trying to learn how to use the check-out and check-in functionality with a file geodatabase.

The database I'm has 13 feature classes in it, most of which are empty because I'm just starting with this project.  The ones that have data are things like the project area, the survey grid, and a transect line.

When I use the 'Get Data for ArcPad' tool on the ArcPad Data Manager toolbar and try to set the action for all the layers I choose 'Checkout all Geodatabase layers and copyout all other layers'.  The problem is that only the line and polygon layers are listed as being checked out:

All the other layers, which are points, are listed as being exported as background data as shapefiles.  All of these feature classes are in the same database, with the same projection.  What is going on with this tool that I can't checkout the other layers?

The option to checkout the point layers isn't even available under action.

Here are the options for the Project Area (polygon) layer:

And here are the options for the Misc_Point layer (point), note that the option for 'Check Out for disconnected editing in ArcPad' is totally missing here:

Any help is appreciated!


texts: Attribute table field default value not transfering to new edited features

texts: Using ESRI ArcMap 10.0, in a versioned environment,  we are in a geo-database and have a feature class with a default value populated in field properties. When we create a new feature using editing tools the default value is not populating the corresponding field as it should.  Has anyone encountered this problem before and any ideas for solutions?


texts: Where can I find a list with beaches which include coordinates?

texts: I use World Weather Online marine API to get wave height of different beaches. This API provide data by coordinates and I need a list with coordinates of beaches.


texts: How to make query of MySQL connection in QGIS?

texts:
i have added a connection with a MySQL database in QGIS
(Add vector layer - Database - MySQL).

Is there an option to make a query over the database in QGIS?


texts: NumPy Array to Feature Class: NameError from Inner Loop Variables

texts:
I have a NumPy Array to FC script that parses a JSON web-service. The outer loop writes to my FC correctly, however my inner loop variables come back as not defined when appending the items.  The outer loop variables include geographic data such as lon/lat, address and other ancillary information about a service request such as service request number, service request type, time created, etc.

The inner loop contains information relating to each service request number, address, etc.

I am using Python 2.6 and ArcGIS 10.2.1

import arcpy
import numpy
import requests
import json
import jsonpickle

fc = "myFC"
f2 =open('C:\Users\GeoffreyWest\Desktop\Request.json')
data2 = jsonpickle.encode(data2)

url2 = "myURL"
headers2 = {'Content-type': 'text/plain', 'Accept': '/'}

r2 = requests.post(url2, data=data2, headers=headers2)
decoded2 = json.loads(r2.text)
print json.dumps(decoded2, sort_keys=True, indent=4)
if arcpy.Exists(fc):
    arcpy.Delete_management(fc)
try:
    r2
except requests.exceptions.ConnectTimeout as e:
    print "Too slow Mojo!"

items = []
for sr in decoded2['Response']['ListOfServiceRequest']['ServiceRequest']:
    SRAddress = sr['SRAddress']
    Latitude = sr['Latitude']
    Longitude = sr['Longitude']
    ReasonCode = sr['ReasonCode']
    SRNumber = sr['SRNumber']
    FirstName = sr['FirstName']
    LastName = sr['LastName']
    ResolutionCode = sr['ResolutionCode']
    HomePhone = sr['HomePhone']
    CreatedDate = sr['CreatedDate']
    UpdatedDate = sr['UpdatedDate']
    BulkyItem = sr['ListOfLa311BulkyItem']
    ElectronicWaste = sr['ListOfLa311ElectronicWaste']
    MoveInMoveOut = sr['ListOfLa311MoveInMoveOut']
    IllegalDumping = sr['ListOfLa311IllegalDumpingPickup']
    ServiceNotComplete = sr['ListOfLa311ServiceNotComplete']
    BrushItems = sr['ListOfLa311BrushItemsPickup']
    Containers = sr['ListOfLa311Containers']
    MHA = sr['ListOfLa311MetalHouseholdAppliancesPickup']
    DeadAnimalRemoval = sr['ListOfLa311DeadAnimalRemoval']
    Manual = sr['ListOfLa311ManualPickup']
    CreatedDate = datetime.datetime.strptime(CreatedDate, "%m/%d/%Y %H:%M:%S")
    UpdatedDate = datetime.datetime.strptime(UpdatedDate, "%m/%d/%Y %H:%M:%S")

    for sr in ElectronicWaste:
            for ewastelocation in ElectronicWaste['La311ElectronicWaste']:
                                    locationewaste =  ewastelocation['CollectionLocation']
    for sr in ElectronicWaste:
            for ewastetype in ElectronicWaste['La311ElectronicWaste']:
                                    itemEwaste =  ewastetype['ElectronicWestType']

 dt = np.dtype([('Address', 'U40'),
        ('LatitudeShape', '<f8'),
        ('LongitudeShape', '<f8'),
        ('Latitude', '<f8'),
        ('Longitude', '<f8'),
        ('ReasonCode','U128'),
        ('SRNumber', 'U40'),
        ('ElectronicWaste', 'U128'),
        ('FirstName', 'U40'),
       ('LastName', 'U40'),
        ('ResolutionCode','U128'),
       ('HomePhone', 'U40'),
        ('CreatedDate', 'U128'),
        ('UpdatedDate', 'U128'),

        ])

    items.append((SRAddress,
                  Latitude,
                 Longitude,
                  Latitude,
                  Longitude,
                  ReasonCode,
                  SRNumber,
                 locationewaste,
                 FirstName,
                  LastName,
                  ResolutionCode,
                  HomePhone,
                  CreatedDate,
                  UpdatedDate,

    ))

sr = arcpy.SpatialReference(4326)

arr = np.array(items,dtype=dt)
NumPyArray = arcpy.da.NumPyArrayToFeatureClass(arr, fc, ['longitudeshape', 'latitudeshape'], sr)
print "success"

Output NameError: name 'locationewaste' is not defined

NOTE: My for loops for electronic waste are indented, I do not know if they are showing up correctly.


texts: Loading features incrementally in a vector layer

texts:
I am new to OpenLayers. In one experiment I am doing I want to load, incrementally, features in a vector layer.
In my setting, these features should be loaded from distinct OSM files (generated dynamically) in a HTTP server.

I know how to create a vector layer from a single OSM file, I have done it with a code similar to this:

...
map = new OpenLayers.Map('map_element', {
    maxResolution: 156543.0399,
    numZoomLevels: 16,
    units: 'm',
    projection: new OpenLayers.Projection('EPSG:900913'),
    displayProjection: new OpenLayers.Projection('EPSG:4326')
});

vectorLayer = new OpenLayers.Layer.Vector("POIs", {
    strategies: [new OpenLayers.Strategy.Fixed()],
    protocol: new OpenLayers.Protocol.HTTP({
    url: "my_dynamically_generated_osm_file.osm",
    format: new OpenLayers.Format.OSM()
    }),
    projection: new OpenLayers.Projection("EPSG:4326")
});

map.addLayers([<other_layers_here>, vectorLayer]);
...

However, I am wondering how I can create first the vector layer without features, and afterwards adding new features defined in different OSM files according to the requests of the user. I know vector layer objects have a method addFeatures(), so I just need to know how could I obtain these features from the distinct OSM files provided by my http server ?


texts: How to convert XYZ to from-to format with azimuth and dip information?

texts:
I have an Excel spreadsheet with header, survey, and assay drill data.

Header data contains hole ID, location coordinates, and length.
Survey data contains related downhole survey with distance, azimuth and dip values.
Assay data contains related location coordinates and %s.

Since I know the coordinates, azimuth, and dip, I would like to be able to convert the assay table to from-to format. Does anyone have a function / procedure / example?

Header data:

Survey data:

Assay data:


texts: Can't find correct coordinate system: UTM83 Zone 15 forcing input in meters; need it to be read in feet

texts:
Forgive me, I'm a CAD person, and am very new to GIS.  I've searched google and this site extensively for an answer, with no luck.

My problem is that I was given a large number of sistered Raster files; a TIFF and an accompanying .tfw file.  The XY coordinates in the .tfw file are given in feet.

I've got these images to come in perfectly in CAD by setting the coordinates to UTM83-15 and manually declaring the input coordinates to be feet.  It lines up perfectly.

However, when I try to use ArcMap, if I set the zone to UTM83 Zone 15N, it seems to force (logically so, I suppose) the input to be read as meters, which pushes everything ~3.28x up and to the right of where I want it.

I cannot for the life of me find a way to force it to be in feet, nor can I find the same coord system not in metric.  Im about to write a script to go through and convert all the X-Y to meteres: an imperfect solution as I don't really want to manipulate my base data.

To fix this:

Is there a way to force ArcMap to read an input file as feet while using a UTM 15 origin?
Is there a foot based surrogate coordinate system with the same origin for UTM83-15N?


texts: How can I add an index to a CartoDB table?

texts:
First of all, do I need to do this or does CartoDB do this based on the query on the fly?

NYC's 311 database is pretty big and it'd be nice to be able to have an index on created_date and complaint_type.


texts:
I need to have quickosm run a query via the modeler.
I have found examples on the net, but I cannot modify them for my own use.
I need to download some vectors and create 3 layers.
How do I do that?
First of all I have to make sure that it always retrieves the data from canvas without indicating the area and then I have to create 3 levels, points, lines and areas.
I am attaching the query code so that I better understand what I need.
qgis 3.16
quickosm 1.14.3
<osm-script output="xml" timeout="25">
    <union>
        <query type="way">
            <has-kv k="railway" v="tram"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="railway" v="rail"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="highway" v="living_street"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="highway" v="pedestrian"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="highway" v="primary"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="highway" v="residential"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="highway" v="secondary"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="highway" v="service"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="highway" v="tertiary"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="highway" v="tertiary_link"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="highway" v="unclassified"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <has-kv k="building"/>
            <bbox-query {{bbox}}/>
        </query>
        <query type="relation">
            <has-kv k="building"/>
            <bbox-query {{bbox}}/>
        </query>
    </union>
    <union>
        <item/>
        <recurse type="down"/>
    </union>
    <print mode="body"/>
</osm-script>
```![enter image description here](https://i.stack.imgur.com/vllsE.jpg)
[![error][1]][1]

  [1]: https://i.stack.imgur.com/gegtP.png


texts: Is there a way to define custom tile names during a gdal_retile.py process, other than the default _01_01. For instance, suppose I want to start naming my tiles at 100_100, is this possible?


texts:
I am trying to upscale a raster using rasterio but I don't know why my process gets killed.

It can print "i am arrived here" but on the reproject() I have as terminal output:

killed

This is the code I am trying to run:

with rasterio.open(mask) as src:
    df_mask = src.read(1)
    aff_mask = src.transform
    profile_mask = src.profile
    width = src.width

profile_esm.update( compress='lzw', dtype = 'float32' )

## using reproject with the same CRS but changing the cell size

new_df= np.empty( shape = ( round(df_mask.shape[0] * 10),
                round(df_mask.shape[1] * 10 )), dtype = 'float32')

newaff = Affine(aff_mask.a / 10, aff_mask.b, aff_mask.c,
                aff_mask.d, aff_mask.e / 10, aff_mask.f)

profile_mask.update( compress='lzw', dtype = 'float32')

print("i am arrived here")
reproject(
        df_mask, new_df,
        src_transform = aff_mask,
        dst_transform = newaff,
        src_crs = src.crs,
        dst_crs = src.crs,
        resampling = Resampling.nearest)

print("i am after")
with rasterio.open('resampling_mask.tif', 'w', **profile_mask) as dst:
    reproject(
        df_mask, new_df,
        src_transform = aff_mask,
        dst_transform = newaff,
        src_crs = src.crs,
        dst_crs = src.crs,
        resampling = Resampling.nearest)

    dst.write(new_df,  indexes=1)


texts: Invalid GeoJSON object, converted from ArcGIS JSON


texts: How to reset the reshape tool in QGIS after 'unknown exception'?

texts:
I am using QGIS 2.18.16 and most of the time I have successfully been using the reshape tool to change the shape of polygons. I have snapping and tracing enabled.

Occasionally when I right-click to finish the reshaping I receive an error 'unknown exception', after which there is no way to undo, reset or start again the tracing. I think the exception is caused because of geometry errors in the polygon that I am reshaping to.

But then I get stuck in reshaping mode. If I close the dialogue window by clicking OK, I am back in the main window and my original reshape trace still exists and I can not find a way to reset it, fix the error, and start again?

All I can do is save the layer, and the project, close down, and reopen QGIS and start again.

Is there an easier and quicker way?


texts: Slow seeding/caching from GeoServer WMS to MapPRoxy

texts:
I face a unormal speed of re-seeding (updating old tiles by new from one source - WMS GeoServer (in GeoServer this data taken from PostgreSQL source)) from 5.32% to this

[08:39:10] 11   5.42% -10175297.20532, 13149614.84996, -10018754.17139, 13306157.88388 (27568 tiles)

reseeding works for about 8 hours!
We use CentOS, 16 gb RAM , SSD.
When seeding was running console throws this information:

  Limited to coverage in: -180.00000, -85.05113, 180.00000, 85.05113 (EPSG:4326)

It is ok? all my data in 3857...
What be the reason of so slowly working (reseeding)? I don't want to update GeoServer to 2.16.0)


texts: Exporting MapSeries using index layer NAME, instead of page number, in ArcPy?

texts:
How can I write this python script to export each page in a MapSeries to a PDF, whilst using the MapSeries index layer NAME as the file name? At the moment, the script is configured to export using the page number.

Is it as simple as replacing 'ms.currentPageNumber' and 'pageNum' with 'ms.currentPageName' and 'pageNam'? And would the {0} need removing on the final line?

import arcpy, os, sys
relpath = os.path.dirname(sys.argv[0])

p = arcpy.mp.ArcGISProject("C:\\Automation\\Automation.aprx")
print ("project located")
l = p.listLayouts()[0]
print ("layouts listed")
if not l.mapSeries is None:
  ms = l.mapSeries
  if ms.enabled:
    for pageNum in range(1, ms.pageCount + 1):
      ms.currentPageNumber = pageNum
      print("Exporting page {0} of {1}".format(str(ms.currentPageNumber), str(ms.pageCount)))
      l.exportToPDF("C:\\Users\\Default\\Production\\Completed_certificates\\Exported_PDF_file_{0}".format(str(ms.currentPageNumber) + ".pdf"))


texts: ArcGIS - Make colour gradient in street network shapefile accordingly to buffer distance from central points (bus stops)

texts:
I need to do a colour gradient in the street network shapefile (cyan selection in the map attached, in a layer) accordingly to the distance from the points in the other layer (the more distant, the darker), always inside the buffer.

In conclusion: I want to get this result (accessibility gradient from bus stops)


texts: Python script for calculating the area of a polygon

texts: I am looking for a python code for this polygon.  I need to calculate the area inside this polygon and I am new to programming. (Very new).  Can anybody assist me with this? I have to use a python code for this, I cannot use any other means.


texts: Multiple Conditional Reclassification - ArcGIS

texts:
This is a question out of curiosity rather than a problem as I am not the most experienced with raster datasets and I wonder if there is a better solution to my current approach. I am trying create a new raster dataset based off the conditions of two raster datasets. i.e. if population size <100 and mean age <10, reclassify to value 1. Now the issue is not creating a conditional statement but if there are several possible categories - say 49, that would require 49 multiple conditional statements. For example:

OutRas = Con((pop_ras < 100) & (age_ras < 10), 1, Con(pop_ras > 100) & ((pop_ras < 200) & (age_ras  < 10), 2, ..., 0))

At the moment I use the following procedure:

1) reclassify population into 7 categories

2) reclassify age into 7 categories

3) add age and population (to create 49 classes)

This approach would require 3 geoprocessing steps as opposed to one but require considerable less time scripting. I feel that there might be a simpler solution that I am not realization
(perhaps a lookup table, SQL). Any thoughts? ideas? Which method is better - speed?


texts:
I am trying to read/write NetCDF file using the NetCDF C++ API:
http://unidata.github.io/netcdf-cxx4/index.html
Part of the program also operates on Shapefile using the GDAL C++ API:
https://gdal.org/tutorials/vector_api_tut.html
I received the following error at line:
NcFile pFile(sFilename_in.c_str(), NcFile::replace);

Warning!  HDF5 library version mismatched error
The HDF5 header files used to compile this application do not match
the version used by the HDF5 library to which this application is linked.
Data corruption or segmentation faults may occur if the application continues.
This can happen when an application was compiled by one version of HDF5 but
linked with a different version of static or shared HDF5 library.
I suspect that the GDAL library uses a different NetCDF library with my NetCDF library.
module load gcc/8.1.0
module load gdal/2.3.1
module load netcdf/4.6.3

Is there any way I check the compatibility of these two? How to resolve this?


texts: How to index USGS ARD Satellite Data into an Open Data Cube?

texts:
I am trying to index a Landsat8-Satellite-Dataset from Massachusetts in the USA. I downloaded the dataset from USGS. Below my question you can find the datacube-metadata.yaml-File. It is surface reflectance data and Analysis Ready Data (ARD). When I am now trying to index the data with: datacube -v dataset add PATH/datacube-metadata.yaml there occurs the error: ERROR No matching Product found for dataset
I think that the problem lies in the Product Definition File. I have used the Product Definition File from Swiss Data Cube: https://github.com/GRIDgva/SwissDataCube/blob/master/ingestors/Landsat/ls8_c1_ard_scene.yaml
I also tried to write one by myself but there is always the same error. Has someone an idea how to solve the error?
datacube-metadata.yaml-File:
acquisition:

aos:
'2018-10-30 15:32:36'

groundstation: {code:
_20}

los:
'2018-10-30 15:33:00'

creation_dt: 2018-10-30 00:00:00

extent:

center_dt:
'2018-10-30 15:32:48'

coord:

ll: {lat:
41.543343249666286,
lon:
-74.26101678236896}

lr: {lat:
41.22740733220001,
lon:
-72.50449089590899}

ul: {lat:
42.84901849362111,
lon:
-73.84246819961363}

ur: {lat:
42.52663821973187,
lon:
-72.05488589831796}

from_dt:
'2018-10-30 15:32:36'

to_dt:
'2018-10-30 15:33:00'

format: {name:
GeoTiff}

grid_spatial:

projection:

geo_ref_points:

ll: {x:
1784415.0,
y: 2264805.0}

lr: {x:
1934415.0,
y: 2264805.0}

ul: {x:
1784415.0,
y: 2414805.0}

ur: {x:
1934415.0,
y: 2414805.0}

spatial_reference:
PROJCS["Albers",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS

84",6378140,298.2569999999957,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]

id: 94675613-a468-45cd-850f-d6423beb4d61

image:

bands:

LINEAGEQA: {path:
LC08_CU_029006_20181030_20190613_C01_V01_LINEAGEQA.tif}

PIXELQA: {path:
LC08_CU_029006_20181030_20190613_C01_V01_PIXELQA.tif}

RADSATQA: {path:
LC08_CU_029006_20181030_20190613_C01_V01_RADSATQA.tif}

SRAEROSOLQA: {path:
LC08_CU_029006_20181030_20190613_C01_V01_SRAEROSOLQA.tif}

SRB1: {path:
LC08_CU_029006_20181030_20190613_C01_V01_SRB1.tif}

SRB2: {path:
LC08_CU_029006_20181030_20190613_C01_V01_SRB2.tif}

SRB3: {path:
LC08_CU_029006_20181030_20190613_C01_V01_SRB3.tif}

SRB4: {path:
LC08_CU_029006_20181030_20190613_C01_V01_SRB4.tif}

SRB5: {path:
LC08_CU_029006_20181030_20190613_C01_V01_SRB5.tif}

SRB6: {path:
LC08_CU_029006_20181030_20190613_C01_V01_SRB6.tif}

SRB7: {path:
LC08_CU_029006_20181030_20190613_C01_V01_SRB7.tif}

satellite_ref_point_end: {x:
29,
y: 6}

satellite_ref_point_start: {x:
29,
y: 6}

instrument: {name:
OLI/TIRS}

lineage:

source_datasets: {}

platform: {code:
LANDSAT_8}

processing_level:
T1

product_type:
LaSRC


texts: Getting "real time" occurence count of values in attribute table for displaying purposes

texts: I'm working with ArcMap 10.2 on a layer of points that I have to move. I'm using the editor to move them. What I also need is to display these points based on how many are sharing the same coordinates X and Y, because that's how I find those I need to move ( I give priority to those points which share coordinates with many other points, for reasons you don't need to be bothered with). Getting a Summary Table from the concatenated X & Y and joining the count back to my points layer is the "standard" way to go, but I need to get a "real time" count of occurence of each value and I cannot extract each time the Summary Table and make the join. I searched for and tried with some python coding in the field calculator but didn't find anything working. Ideas?


texts: Clipping soil survey map service to feature class?

texts:
I'm guessing this is not possible but is there any way to clip a Esri map service such as the Soil Survey Map for a another feature polygon? http://server.arcgisonline.com/ArcGIS/rest/services/Specialty/Soil_Survey_Map/MapServer

I know it's a not a feature service but what is confusing is how if you zoom in far enough it will show the soil series information.  What I would like to do is clip it to a smaller area like less than 1000 acres based on a polygon feature class.  Maybe someone knows of another feature service available that would show the soil survey down to the soil series detail like this Esri map service does?


texts:
I have a feature class in which one field (let's call it "Field_A"), as yet empty, is to be populated by the contents of another field ("Field_B") multiplied by a certain number (I'll call it the "constant"). Would it be a good idea to simply populate Field_A using an update cursor, as follows:

rows = arcpy.UpdateCursor(InTable)
for row in rows:
     row.Field_A = row.Field_B * constant
     rows.updateRow(row)
del row, rows

, or do I have to use the field calculator? If I have to use the field calculator, I'm having an issue with the syntax. I defined the expression thus:

expression = "Field_A * constant"

and the code is as follows:

arcpy.CalculateField_management(InTable, "Field_A", expression)

When using CalculateField as shown, the script runs without errors, but "Field_A" is populated with zeros as a result. What am I doing wrong with CalculateField, or should I simply use an update cursor?


texts: Finding which census tracts lie within a given shape?

texts:
I'm new to GIS in R, and I am trying to create a crosswalk that maps the communities drawn in these shapefiles to census tract shapefiles for Los Angeles County that I've gathered from the Tigris package.
What I'd like is a list of census tracts that correspond to each geography--for instance, Agoura Hills contains x,y,z tracts, Alhambra contains a,b,c tracts, etc. The issue is that some of these places aren't cities but neighborhoods for larger cities which is why I can't just use total census to get what I need.
I'm not entirely sure, but I think what I might be trying to do is a spatial join?


texts: Creating a 3D model of airspace that I can interrogate to see if some other feature (like a tree) will intrude into the airspace

texts:
The goal is to create a 3D model of imaginary planes in space that I can interrogate to see where other objects in the vicinity of the model cut through any part of that plane. For example, an airspace surrounding an airport and a tall tree that extends into the lower surface of the airspace.

The model will consist various horizontal planes and sloped planes that exist between the horizontal planes (kind of like a dinner plate shape).

I thought I could try creating a DEM that I could analyse. I tried to make polygons around a point (buffer), then converted them to lines, then to points - because the interpolation (average) tool required scattered data (i.e. points) with z values). It didn't work - the output raster wasn't to scale. I also tried using the Interpolation Plugin and the output was to scale but it's not precise enough for the purpose. The dimensions of the airspace are mathematically defined using gradients, distances and elevations.

I created polygons easily enough with the Z values appended to it that represent horizontal surfaces. I just can't figure out how to create the sloped surfaces with varying Z values (which depend on the distance from some point in space, like the edge of a horizontal surface polygon or a point).

Z values need to be based on height above an earth reference point, such as metres above Australian Height Datum (AHD). So, I'm not sure if Sketchup will be useful here (???). I need to be able to interrogate the model to see if an object in Earth is going to penetrate the surface.

Does anyone have any suggestions regarding this issue - or generally some other way to go about this problem using another program?


texts:
I was given a json file that contains information like as below (truncated version but complete in structure)-

{"type": "FeatureCollection", "features": [{"geometry": {"type": "Polygon", "coordinates": [[[3909262.46926247, -1770881.94706788], [3909261.06123541, -1770740.38745684], [3909401.80402645, -1770738.96964834], [3909400.38998219, -1770596.88864125], [3909400.38998219, -1770596.88864125], [3909539.73283601, -1770453.91231545], [3909539.73283601, -1770453.91231545], [3909824.04203656, -1770734.7109451], [3909404.6218828, -1771022.08945481], [3909263.87730952, -1771023.50751348], [3909262.46926247, -1770881.94706788]]]}, "type": "Feature", "properties": {"dtegranted": 1321488000000, "iss_date": "17/11/2011", "objectid": 1, "company": "Terrastone Limited", "lic_no": "ML0002", "layer": "Mining Licence", "exp_date": "17/11/2016", "address": "P.O. BOX 995, LILONGWE", "id": 1, "area_sqkm": 0.169516104612393, "district": "NJULI", "license": "ML0002", "area": "169516.104612393", "minerals": "Rock Aggregate", "contact": "CHILADZULU", "dteexpires": 1479340800000, "type": "Mining Licence", "statusgrp": "Active"}}]}

and said that projection format is UTM36S nothing else. Now I need to convert this data back into WGS84 and create a polygon shapefile.

I am failing to convert into WGS84 and I converted those json files into shapefile using global mapper(used UTM36S during export to shapefile) since arcmap fails but that pops up as below though global mapper was able decode attributes correctly from the json file but messes up with the geometry-

I tried to convert those utm x and y data into WGS84 using python function as mentioned here but it gives me wrong coordinate
Please help-

Full version of data is at http://pastebin.com/embed_js/9FfDYxgN.


texts: Using R with ArcGIS for Desktop (without using GME)?

texts:
I am trying to learn R generally and making statistics for maps. I found some useful tricks for making maps with it but it is not the best you can do. My everyday GIS program is ArcGIS.  How I can use R in ArcGIS?  A tool I use is the GME plugin. Is this the only way?

I tried to use R in QGIS but a plugin I found never worked.